(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* ::Input::Initialization:: *)
(* Start package *)
BeginPackage["MaXrd`"];

(* Import usage messages from file *)
If[!FailureQ@FindFile[#],Get@FindFile[#]]&["MaXrd/Core/Messages.txt"];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
AttenuationCoefficient::invalidSource="The source \[LeftGuillemet]`1`\[RightGuillemet] is invalid.";
AttenuationCoefficient::invcoeff="The \[LeftGuillemet]`1`\[RightGuillemet] coefficient is not recognised.";
AttenuationCoefficient::peOnly="Only the photoelectric cross section is related to \!\(\*
StyleBox[FormBox[SuperscriptBox[\"f\", \"\[Prime]\[Prime]\",\nMultilineFunction->None],
TraditionalForm], \"TI\"]\).";
AttenuationCoefficient::asfLambda="The wavelength, `1` \[CapitalARing], must be smaller than 2.5 \[CapitalARing] when using \!\(\*FormBox[SuperscriptBox[\(f\), \(\[DoublePrime]\)],
TraditionalForm]\).";

Options@AttenuationCoefficient={
"Coefficient"->"LinearAttenuation",
"MassCoefficientMethod"->"DivideByDensity",
(* GetScatteringCrossSections *)
"PhysicalProcess"->"",
"Source"->"xraylib",
"Units"->True
};

SyntaxInformation@AttenuationCoefficient={
"ArgumentsPattern"->{_,_.,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
AttenuationCoefficient[
crystal_String,
lambda:_?(NumericQ[#]||QuantityQ[#]&):-1,
OptionsPattern[]]:=Block[{
\[Lambda],unitsQ=TrueQ@OptionValue["Units"],
coeff=OptionValue["Coefficient"],
mcMethod=OptionValue["MassCoefficientMethod"],
source=OptionValue["Source"],pp=OptionValue["PhysicalProcess"],csDir,asfFile,
V,\[Sigma],\[CapitalSigma],column,atomdata,r,siteM,
fpp,re,\[Mu],\[Rho]
},

(*---* Checking input *---*)
	InputCheck["CrystalQ",crystal];
	\[Lambda]=InputCheck["ProcessWavelength",crystal,lambda];

	csDir=FileNameJoin[{
	$MaXrdPath,"Core","Data","CrossSections",source}];
	asfFile=FileNameJoin[{$MaXrdPath,"Core","Data",
	"AtomicScatteringFactor","AnomalousCorrections",source<>".m"}];
	If[DirectoryQ[csDir]\[Nor]FileExistsQ[asfFile],
		Message[AttenuationCoefficient::invalidSource,source];
		Abort[]];

	(* Processing cross section type *)
	If[
	(* a. Select 'PhysicalProcess' manually *)
	pp=!="",
	column=Which[
	MemberQ[{
	"Photoelectric","Photoionisation"},pp],2,
	MemberQ[{
	"Coherent","Rayleigh","Thompson",
	"Classical","Elastic"},pp],3,
	MemberQ[{
	"Incoherent","Compton","Inelastic"
	},pp],4,
	pp==="Total",5,
	True,Message[AttenuationCoefficient::invproc,pp];Abort[]
	],

	(* b. Select automatically based on coefficient type *)
	column=Which[
	MemberQ[
	{"LinearAttenuation","MassAbsorption"},
	coeff],
	5,(* Total = ph.el. + Ray. + Comp. *)

	True,
	Message[AttenuationCoefficient::invcoeff,coeff];
	Abort[]]
	];

(*---* Calculations *---*)
	(* Auxiliary variables *)
	V=Sqrt@Det@GetCrystalMetric@crystal;
	
	(* Calculation method *)
	If[DirectoryQ@csDir,
	(* a. Using cross sections *)
	\[Sigma]=Normal@GetScatteringCrossSections[crystal,\[Lambda],
	"PhysicalProcess"->pp,"Source"->source,"Units"->False];

		(* Multiplying atoms with corresponding cross sections *)
		atomdata=Values@$CrystalData[[crystal,"AtomData",All,
	{"Element","FractionalCoordinates","OccupationFactor"}]];
		atomdata[[All,1]]=StringDelete[atomdata[[All,1]],
		{"+","-",DigitCharacter}];
		atomdata=atomdata/.Join[
		\[Sigma],{Missing["KeyAbsent","OccupationFactor"]->1}];
		r=atomdata[[All,2]];
		siteM=Table[Length@SymmetryEquivalentPositions[
		crystal,r[[a]],"UseCentring"->True],{a,Length@r}];
		atomdata[[All,2]]=siteM;
		\[CapitalSigma]=Total[Times@@@atomdata];
		\[Mu]=\[CapitalSigma]/V,

	(* b. Using f-double-prime *)
		(* Check wavelength *)
		If[\[Lambda]>2.5,
		Message[AttenuationCoefficient::asfLambda,ToString@\[Lambda]];
		Abort[]];

		column=3;(* Force p.e. cross section only *)
		re=2.81794032*^-15;
		fpp=Im@StructureFactor[crystal,{0,0,0},\[Lambda],
		"AbsoluteValue"->False,
		"f1f2Source"->source];
		(* (See formula in documentation page details) *)
		\[Mu]=2*re*\[Lambda]/V*Abs[fpp]*Power[10,18]
	];

	(* Normalise by mass density? *)
	If[coeff==="MassAbsorption"&&mcMethod==="DivideByDensity",
	\[Rho]=CrystalDensity[crystal,"Units"->False];
	If[unitsQ,
	Return@Quantity[\[Mu]/\[Rho],"Centimeters"^2/"Grams"],
	Return[\[Mu]/\[Rho]]]
	];

(* Output linear coefficient *)
If[unitsQ,Quantity[\[Mu],"Centimeters"^(-1)],\[Mu]]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
BraggAngle::invinput="Input must either be the name of a crystal or a metric matrix.";

Options@BraggAngle={
"Units"->True,
"AngleThreshold"->90.*Degree
};

SyntaxInformation@BraggAngle={
"ArgumentsPattern"->{_,_.,_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
BraggAngle[
input_,
lambda:_?(NumericQ[#]||QuantityQ[#]&):-1,
reflections_List,
OptionsPattern[]]:=Block[{
hkl,G,H,\[Lambda]=N@lambda,sl,bragg,angle,angleThreshold
},
(* Check crystal (or metric) and reflection(s) *)
Which[
StringQ@input,
	InputCheck["CrystalQ",input];
	G=GetCrystalMetric[input],
MatrixQ@input,
	G=input,
True,
	Message[BraggAngle::invinput];Abort[]
];

hkl=InputCheck[reflections,"Integer","WrapSingle"];

(* Reciprocal metric *)
H=Chop@N@Inverse@G;

(* Process wavelength *)
\[Lambda]=If[MatrixQ@input,
InputCheck["GetEnergyWavelength",\[Lambda],False],
InputCheck["ProcessWavelength",input,\[Lambda]]
];

(* Sin/lambda, from Bragg's law and inner product *)
sl[h_]:=Sqrt[h.H.h]/2.;
bragg[h_]:=N[ArcSin[sl[h]*\[Lambda]]]/.x_Complex->Undefined;

(* Bragg's law *)
angle=bragg/@hkl;

(* Optional: Truncate at chosen angle threshold *)
angleThreshold=OptionValue["AngleThreshold"];
If[0<=angleThreshold<\[Pi]/2,
angle=Select[angle,(#<=angleThreshold)&]
];
angle=angle/Degree;

(* Option: Units *)
If[OptionValue["Units"],
angle=Quantity[angle,"Degrees"];
angle=angle/.Quantity[Undefined,"Degrees"]->Undefined
]; 

(* If only one reflection, return content *)
If[Length@angle==1,
First@angle,angle]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
ConstructDomains::SectorRegionsInvalidNumberOfPairs="The number of pairs must be a natural number.";
ConstructDomains::SectorRegionsInvalidWidth="Angular width of the sectors must be a number";

Options@ConstructDomains={
"TransitionProbabilities"-><|
0->0.95,1->0.92,2->0.86,3->0.75,4->0.40,5->0.50,6->0.75,7->0.12,8->0.03|>
};

SyntaxInformation@ConstructDomains={
"ArgumentsPattern"->{_,_,_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
ConstructDomains[
{A_Integer,B_Integer,C_Integer},
numberOfDomains_Integer,numberOfCycles_Integer,
OptionsPattern[]]:=Block[{
structureSize=A*B*C,
insertionCoordinates,domainTable,transitionProbabilities=OptionValue["TransitionProbabilities"],
visitOrderInCurrentCycle,currentCellIndex,currentCell,nearest,nearestFiltered,neighbourDomains,numberOfDomainsEqualToSelf
},

insertionCoordinates=Flatten[Table[{i,j,k},
{i,0,A-1},{j,0,B-1},{k,0,C-1}],2];
domainTable=Association@Thread[insertionCoordinates->RandomSample[
Flatten@ConstantArray[
Range@numberOfDomains,\[LeftCeiling]structureSize/numberOfDomains\[RightCeiling]],
structureSize]];

Do[
visitOrderInCurrentCycle=RandomSample@Range@structureSize;
Do[
currentCellIndex=visitOrderInCurrentCycle[[i]];
currentCell=insertionCoordinates[[currentCellIndex]];
nearest=Nearest[insertionCoordinates,currentCell,26];
nearestFiltered=Select[
Delete[nearest,1],
SquaredEuclideanDistance[currentCell,#]<=3.&];
neighbourDomains=Lookup[domainTable,nearestFiltered,0.];
numberOfDomainsEqualToSelf=Count[neighbourDomains,domainTable[currentCell]];
If[Random[]<transitionProbabilities[numberOfDomainsEqualToSelf],
domainTable[currentCell]=RandomChoice@neighbourDomains],
{i,3}],numberOfCycles
];

(* Express domain table as '{outputSize, {domains}}' *)
{{A,B,C},Values@domainTable}
]


(* ::Input::Initialization:: *)
ConstructDomains["SectorRegions",
{A_Integer,B_Integer,C_Integer},regionSettings_List,
OptionsPattern[]]:=Block[{
insertionCoordinates,identifiers,
regions,
MakeSectorRegion,angleStarts,oppositeStarts,allStarts,angleRanges,disks,
FindMatch,n
},

(* Input cheks *)
If[!AllTrue[regionSettings[[All,1]],IntegerQ[#]&&Positive[#]&],
Message[ConstructDomains::SectorRegionsInvalidNumberOfPairs];Abort[]];
If[!AllTrue[regionSettings[[All,2]],NumericQ],
Message[ConstructDomains::SectorRegionsInvalidNumberOfPairs];Abort[]];

(* Auxiliary functions *)
MakeSectorRegion[
numberOfPairs_Integer,
width_?NumericQ,
startAngle_:Null]:=(

If[NumericQ@startAngle,
angleStarts=Table[a,{a,startAngle,2.*\[Pi],2.*\[Pi]/(numberOfPairs*2)}],
angleStarts=RandomReal[{0.,2.*\[Pi]},numberOfPairs]
];
oppositeStarts=angleStarts+\[Pi];
allStarts=Join[angleStarts,oppositeStarts];
angleRanges={#,#+width}&/@allStarts;
disks=Disk[{A/2,B/2},A,#]&/@angleRanges;
RegionUnion@@disks
);

FindMatch[{x_,y_,z_}]:=(
n=Do[If[RegionMember[regions[[i]],{x,y}],Return@i],
{i,Length@regions}];
{x,y,z}->n);

(* Main procedure *)
insertionCoordinates=Flatten[Table[{i,j,0},
{i,0,A-1},{j,0,B-1}],1];
regions=MakeSectorRegion@@@regionSettings;
insertionCoordinates=(FindMatch/@insertionCoordinates)
/.{Null->Length@regions+1};
identifiers=Values@insertionCoordinates;

If[C>1,
identifiers=Flatten@Transpose@ConstantArray[identifiers,C]
];

(* Express domain table as '{outputSize, {domains}}' *)
{{A,B,C},identifiers}
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
Options@CrystalDensity={
"Units"->True
};

SyntaxInformation@CrystalDensity={
"ArgumentsPattern"->{_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
CrystalDensity[crystal_String,OptionsPattern[]]:=Block[
{data,unitsQ,Z,
f,m,V,NA,
X,o,xyz,M,mass,
temp},

(*---* Input check *---*)
	InputCheck["CrystalQ",crystal];
	data=$CrystalData[crystal];
	unitsQ=OptionValue["Units"];

	(* Return density if contained in '$CrystalData' *)
	temp=data["MassDensity"];
	If[!MissingQ[temp],Return@temp];

	(* Lookup formula units *)
	Z=Lookup[data,"FormulaUnits",-1];

(*---* Calculations *---*)
(*--* Common variables *--*)
	(* Volume, converted to cm^3 *)
	V=Sqrt@Det@GetCrystalMetric@crystal;
	If[unitsQ,
	V=Quantity[V*10^(-24),"Centimeters"^3],
	V=V*10^(-24)];

	(* Chemical formula *)
	f=Sort@GetElements[crystal,"Tally"->True];

If[Positive@Z,
(*--* A. Calculate \[Rho] from Z *--*)
	(* Atomic mass of one unit *)
	m=MapAt[$PeriodicTable[#,"AtomicMass"]&,f,{All,1}];
	m=Total[Times@@#&/@m];
	If[unitsQ,m=Quantity[m,"Grams"/"Moles"]];

	(* Avogadro's constant *)
	If[unitsQ,
	NA=UnitConvert[Quantity["AvogadroConstant"],"Moles"^-1],
	NA=6.022140857*^23];

	(* Calculating \[Rho] *)
	Return[(Z*m)/(V*NA)],


(*--* B. Calculate \[Rho] from atom data, symmetry and occupation *--*)
	(* Elements, occupation factors and coordinates *)
	{X,o,xyz}=Transpose@Values@data[["AtomData",All,
{"Element","OccupationFactor","FractionalCoordinates"}]];
	X=StringDelete[X,{DigitCharacter,"+","-"}];
	o=o/._Missing->1.;

	(* Site multiplicities *)
	M=o*(Length/@(SymmetryEquivalentPositions[
data["SpaceGroup"],xyz]));

	(* Counting *)
	temp=Transpose/@GatherBy[Transpose[{X,M}],First];
	temp=Sort[temp/.{x_List,m_List}/;
Depth[x]===2:>{First@x,Total@m}];

	(* Total atom mass *)
	mass=Total[temp/.{X_String,f_?NumericQ}:>
	f*$PeriodicTable[X,"AtomicMass"]];
	If[unitsQ,
	mass=UnitConvert[Quantity[mass,"AtomicMassUnit"],"Grams"],
	mass=mass*(1.6605390*^-24)];

	(* Calculated density *)
	Return[mass/V]
]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
CrystalFormulaUnits::mismatch="Element mismatch detected.";

Options@CrystalFormulaUnits={
"IgnoreHydrogen"->True
};

SyntaxInformation@CrystalFormulaUnits={
"ArgumentsPattern"->{_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
CrystalFormulaUnits[
crystal_String,OptionsPattern[]]:=Block[
{data,
\[Rho],f,fWithoutH,m,V,NA,
X,XwithoutH,o,xyz,M,Z,
temp},

(*---* Input check *---*)
	InputCheck["CrystalQ",crystal];
	data=$CrystalData[crystal];

	(* Return Z if contained in '$CrystalData' *)
	temp=data["FormulaUnits"];
	If[!MissingQ[temp],Return@temp];

	(* Lookup density *)
	\[Rho]=Lookup[data,"MassDensity",-1];

(*---* Calculations *---*)
	(* Chemical formula *)
	f=Sort@GetElements[crystal,"Tally"->True];

	(* Atomic mass of one unit *)
	m=MapAt[$PeriodicTable[#,"AtomicMass"]&,f,{All,1}];
	m=Total[Times@@#&/@m];

	(* Volume [cubic centimeters] and Avogadro's constant *)
	V=Sqrt@Det@GetCrystalMetric[crystal]*Power[10,-24];
	NA=6.022140857*^23;

If[Positive@\[Rho],
(*--* A. Calculate Z from mass denisty *--*)
	(* Calculate formula units from density *)
	If[QuantityQ@\[Rho],\[Rho]=QuantityMagnitude@
	UnitConvert[\[Rho],"Grams"/"Centimeters"^3]];

	(* Calculate formula units from density *)
	Return[\[Rho]*V*NA/m],


(*--* B. Calculate Z by counting symmetry-generated atoms *--*)
	(* Elements, occupation factors and coordinates *)
	{X,o,xyz}=Transpose@Values@data[["AtomData",All,
{"Element","OccupationFactor","FractionalCoordinates"}]];
	X=StringDelete[X,{DigitCharacter,"+","-"}];
	XwithoutH=DeleteCases[X,"H"];
	o=o/._Missing->1.;
	fWithoutH=Sort@DeleteCases[f[[All,1]],"H"];
	
		(* Check: Compare formula and atom data *)
		If[!TrueQ@OptionValue["IgnoreHydrogen"],
		fWithoutH=f;XwithoutH=X];
		If[Sort@DeleteDuplicates@XwithoutH=!=fWithoutH,
		Message[CrystalFormulaUnits::mismatch];Abort[]];

	(* Site multiplicities *)
	M=o*(Length/@(SymmetryEquivalentPositions[
data["SpaceGroup"],#]&/@xyz));

	(* Counting *)
	temp=Transpose/@GatherBy[Transpose[{X,M}],First];
	temp=Sort[temp/.{x_List,m_List}/;
Depth[x]===2:>{First@x,Round@Total@m}];

	(* Check if hydrogen is ignored *)
	If[MemberQ[f,{"H",_}]&&(!MemberQ[temp,{"H",_}]),
	f=DeleteCases[f,{"H",_}];
	temp=DeleteCases[temp,{"H",_}]];

	(* Check agreement of 'Z' *) 
	Z=temp[[All,2]]/f[[All,2]];
		(* a. Common integer factor *)
		If[MatchQ[DeleteDuplicates@Z,{_Integer}],
		Return@First@Z,

		(* b. Calculate density, then find Z *)
		\[Rho]=CrystalDensity[crystal,"Units"->False];
		Z=(\[Rho]*V*NA)/(m);
		Return@Z
	]
]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
CrystalPlot::InvalidMethod="\"AtomSizeMethod\" must either be \"RadiusTable\" or \"DisplacementParameters\".";
CrystalPlot::InvalidDisplay="\"UnitCell\" must be either \"Box\", \"Axes\" or \"None\".";

Options@CrystalPlot=SortBy[Normal@Union[
Association@Options@Graphics3D,<|
"AxisFunction"->Line,
"OpacityMap"-><||>,
"RGBStyle"->True,
"StructureSize"->{0,0,0},
"UnitCellDisplay"->"Box",
(* Graphics3D *)
Boxed->False,
PlotRange->All,
Lighting->"Neutral",
SphericalRegion->True
|>],ToString[#[[1]]]&];

SyntaxInformation@CrystalPlot={
"ArgumentsPattern"->{_,OptionsPattern[{CrystalPlot,Graphics3D}]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
CrystalPlot[
crystalInput_String,
OptionsPattern[{CrystalPlot,Graphics3D}]]:=Block[{
crystalDataOriginal=$CrystalData,
crystal=crystalInput,structureSize=OptionValue["StructureSize"],
rgbStyle=TrueQ@OptionValue["RGBStyle"],latticeStyleList,CreateBoxEdges,toCartesianMatrix,MakeSpheres,
latticePlotFunction=OptionValue["AxisFunction"],
crystalData,atomData,coordinates,
atomDataTemp,spheres,
basisArrowCoordinates,translations,coordinatePairs,unitCellPlotData,unitCellDisplay=OptionValue["UnitCellDisplay"],
opacityMap=OptionValue["OpacityMap"],
plotOptions
},

(* Optional: Expand crystal *)
If[structureSize=!={0,0,0},
ExpandCrystal[crystalInput,structureSize,"StoreTemporarily"->True];
crystal=(Keys@MaXrd`Private`$TempCrystalData)[[1]];
$CrystalData=MaXrd`Private`$TempCrystalData;
];

(* Auxiliary *)
If[!ValueQ@MaXrd`Private`atomRadiusTable,
MaXrd`Private`atomRadiusTable=<|"H"->0.53,"He"->0.31,"Li"->1.67,"Be"->1.12,"B"->0.87,"C"->0.67,"N"->0.56,"O"->0.48,"F"->0.42,"Ne"->0.38,"Na"->1.9,"Mg"->1.45,"Al"->1.18,"Si"->1.11,"P"->0.98,"S"->0.87,"Cl"->0.79,"Ar"->0.71,"K"->2.43,"Ca"->1.94,"Sc"->1.84,"Ti"->1.76,"V"->1.71,"Cr"->1.66,"Mn"->1.61,"Fe"->1.56,"Co"->1.52,"Ni"->1.49,"Cu"->1.45,"Zn"->1.42,"Ga"->1.36,"Ge"->1.25,"As"->1.14,"Se"->1.03,"Br"->0.94,"Kr"->0.87,"Rb"->2.65,"Sr"->2.19,"Y"->2.12,"Zr"->2.06,"Nb"->1.98,"Mo"->1.9,"Tc"->1.83,"Ru"->1.78,"Rh"->1.73,"Pd"->1.69,"Ag"->1.65,"Cd"->1.61,"In"->1.56,"Sn"->1.45,"Sb"->1.33,"Te"->1.23,"I"->1.15,"Xe"->1.08,"Cs"->2.98,"Ba"->2.53,"La"->1.5,"Ce"->1.5,"Pr"->2.47,"Nd"->2.06,"Pm"->2.05,"Sm"->2.38,"Eu"->2.31,"Gd"->2.33,"Tb"->2.25,"Dy"->2.28,"Ho"->2.26,"Er"->2.26,"Tm"->2.22,"Yb"->2.22,"Lu"->2.17,"Hf"->2.08,"Ta"->2.,"W"->1.93,"Re"->1.88,"Os"->1.85,"Ir"->1.8,"Pt"->1.77,"Au"->1.74,"Hg"->1.71,"Tl"->1.56,"Pb"->1.54,"Bi"->1.43,"Po"->1.35,"At"->1.27,"Rn"->1.2,"Fr"->1.5,"Ra"->1.5,"Ac"->1.5,"Th"->1.5,"Pa"->1.5,"U"->1.5,"Np"->1.5,"Pu"->1.5,"Am"->1.5,"Cm"->1.5,"Bk"->1.5,"Cf"->1.5,"Es"->1.5,"Fm"->1.5,"Md"->1.5,"No"->1.5,"Lr"->1.5,"Rf"->1.5,"Db"->1.5,"Sg"->1.5,"Bh"->1.5,"Hs"->1.5,"Mt"->1.5,"Ds"->1.5,"Rg"->1.5,"Cn"->1.5,"Nh"->1.5,"Fl"->1.5,"Mc"->1.5,"Lv"->1.5,"Ts"->1.5,"Og"->1.5|>];

latticeStyleList=ConstantArray[Black,12];
If[rgbStyle,latticeStyleList[[;;3]]={Red,Green,Blue}
];

CreateBoxEdges[{a_,b_,c_},{t1_,t2_,t3_}]:={
a,b,c,
t2[a],t1[b],
t1[c],t2[c],t1[t2[c]],
t3[a],t3[b],t3[t2[a]],t3[t1[b]]
};

toCartesianMatrix=Transpose@GetCrystalMetric[
crystal,"ToCartesian"->True];

MakeSpheres[element_,coordinateList_]:={
ColorData["Atoms"][element],
If[opacityMap=!=<||>,Opacity[Lookup[opacityMap,element,1.0]]],
Sphere[
Transpose[toCartesianMatrix].#&/@coordinateList,
MaXrd`Private`atomRadiusTable@element
]
};

(* Acquiring data *)
InputCheck["CrystalQ",crystal];
crystalData=$CrystalData[crystal];
atomData=crystalData["AtomData"];
atomData[[All,"Element"]]=StringDelete[
atomData[[All,"Element"]],
{"+","-",DigitCharacter}];
coordinates=atomData[[All,"FractionalCoordinates"]];

(* Preparing atom spheres *)
atomDataTemp=GatherBy[Lookup[atomData,
{"Element","FractionalCoordinates"}],#[[1]]&];
atomDataTemp=AssociationThread[
atomDataTemp[[All,1,1]]->atomDataTemp[[All,All,2]]];
spheres=KeyValueMap[MakeSpheres,atomDataTemp];

(* Basis/lattice vectors and boxes *)
basisArrowCoordinates={{0,0,0},#}&/@toCartesianMatrix;
translations=TranslationTransform/@toCartesianMatrix;

Which[
unitCellDisplay==="Box",
coordinatePairs=CreateBoxEdges[
basisArrowCoordinates,translations];
unitCellPlotData=Table[{
latticeStyleList[[i]],
If[i>3,#,Arrow@#]&@latticePlotFunction@coordinatePairs[[i]]
},{i,Length@coordinatePairs}],

unitCellDisplay==="Axes",
unitCellPlotData=Transpose[{
If[rgbStyle,{Red,Green,Blue},ConstantArray[Black,3]],
Arrow[latticePlotFunction[{{0,0,0},#}]]&/@toCartesianMatrix
}],

unitCellDisplay==="None",unitCellPlotData={},

True,Message[CrystalPlot::InvalidDisplay];Abort[]
];

(* If crystal was expanded, reset pointer to original $CrystalData *)
If[structureSize=!={0,0,0},
$CrystalData=crystalDataOriginal
];

(* Plot options *)
plotOptions=Association@FilterRules[
#->OptionValue[#]&/@(Keys@Options@CrystalPlot),
Options@Graphics3D];

If[
MemberQ[{"Trigonal","Hexagonal"},
GetSymmetryData[crystal,"CrystalSystem"]]&&OptionValue["ViewPoint"]===OptionValue[Graphics3D,ViewPoint],
AssociateTo[plotOptions,ViewPoint->Above]
];

(* Plot *)
Graphics3D[
Join[unitCellPlotData,spheres],
Sequence@@Normal@plotOptions]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
Options@DarwinWidth={
"Units"->True,
	(* ExtinctionLength *)
	"Polarisation"->"\[Pi]"
};

SyntaxInformation@DarwinWidth={
"ArgumentsPattern"->{_,_,_.,_.,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
DarwinWidth[
crystal_,
lambda:_?(NumericQ[#]||QuantityQ[#]&):-1,
hklInput_List,
\[Gamma]o:_?NumericQ:1.0,
\[Gamma]h:_?NumericQ:1.0,
OptionsPattern[]]:=Block[{
hkl,L,\[Lambda],\[Theta],\[CapitalLambda]o,\[Delta]os},

(* Check input *)
InputCheck["CrystalQ",crystal];
hkl=InputCheck[hklInput,"Integer","WrapSingle"];
L=Length@hkl;
\[Lambda]=InputCheck["ProcessWavelength",crystal,lambda];

(* Miscellaneous preparations *)
\[Theta]=BraggAngle[crystal,\[Lambda],hkl,"Units"->False]*Degree;
\[CapitalLambda]o=ExtinctionLength[crystal,\[Lambda],hkl,\[Gamma]o,\[Gamma]h,
"Polarisation"->OptionValue["Polarisation"],
"Units"->False];

(* Darwin width *)
\[Delta]os=2*\[Lambda]/\[CapitalLambda]o*\[Gamma]h/Sin[2\[Theta]]*(1 (* \[Mu]m *)/(10^4) (* \[CapitalARing] *))*((10^6) (* \[Mu]rad *)/1 (* rad *));
\[Delta]os=Flatten[{Chop@\[Delta]os}];

(* Option: Units *)
If[OptionValue["Units"],
\[Delta]os=\[Delta]os/.x_?NumericQ:>Quantity[x,"Microradians"]];

If[L===1,First@\[Delta]os,\[Delta]os]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
DistortStructure::InvalidDistortionType="\"DistortionType\" must be set to either \"Crystallographic\" or \"Cartesian\".";
DistortStructure::InvalidDimensions="Function and variables must both be three-dimensional.";

Options@DistortStructure={
"DistortionType"->"Crystallographic",
"NewLabel"->""
};

SyntaxInformation@DistortStructure={
"ArgumentsPattern"->{_,_,_,OptionsPattern[{CrystalPlot,VectorPlot3D}]},
"LocalVariables"->{"Solve",{3}}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
DistortStructure[crystal_,vectorField_,variables_,OptionsPattern[]]:=Block[{
newLabel=OptionValue["NewLabel"],
distortionType=OptionValue["DistortionType"],
M,inverseM,f,distortions,
coordinates,coordinatesCartesian,newCoordinates,
crystalCopy
},

(* Input checks *)
InputCheck["CrystalQ",crystal];
If[!StringQ@newLabel||newLabel==="",newLabel=crystal];
If[!MemberQ[{"Cartesian","Crystallographic"},distortionType],
Message[EmbedStructure::InvalidDistortionType];Abort[]];
If[Length/@{vectorField,variables}=!={3,3},
Message[DistortStructure::InvalidDimensions];Abort[]];

(* Vector field function *)
f[xyz_]:=N@vectorField/.Thread[variables->xyz];

(* Calculate individual distortions *)
coordinates=$CrystalData[[crystal,"AtomData",All,"FractionalCoordinates"]];
distortions=f/@coordinates;

(* Determine distortion type *)
If[distortionType==="Cartesian",
M=GetCrystalMetric[crystal,"ToCartesian"->True];
inverseM=Inverse@M;
coordinatesCartesian=M.#&/@coordinates;
newCoordinates=coordinatesCartesian+distortions;
newCoordinates=inverseM.#&/@newCoordinates,

(* Shifts in a pure crystallographic frame *)
newCoordinates=coordinates+distortions
];

(* Create new entry in `$CrystalData` *)
crystalCopy=$CrystalData[crystal];
crystalCopy[["AtomData",All,"FractionalCoordinates"]]=newCoordinates;
AssociateTo[$CrystalData,newLabel->crystalCopy];

(* Update auto-completion *)
InputCheck["Update$CrystalDataAutoCompletion"];

newLabel
];


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
DomainPlot::InputMismatch="Input size does not match size of domain list.";
DomainPlot::InvalidDomainIndex="Domain representations must be non-negative integers";

Options@DomainPlot={
"Colours"->{Red,Green,Blue,Yellow,Purple},
"CrystalFamily"->"Cubic",
"GraphicFunction"->Automatic,
Opacity->1.0,
"RotationAnchorReference"->"DomainCentroid",
"RotationAnchorShift"->{0,0,0},
"RotationMap"-><||>
};

SyntaxInformation@DomainPlot={
"ArgumentsPattern"->{{{_,_,_},_},OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
DomainPlot[
{{A_Integer,B_Integer,C_Integer},domains_List},
OptionsPattern[]]:=Block[{
twoDimensionalQ,crystalFamily=OptionValue["CrystalFamily"],graphicFunction=OptionValue["GraphicFunction"],coordinates,coordinateDomainMap,M,makePolytope,
rotationMap=OptionValue["RotationMap"],rotationQ,R,
anchorShift=OptionValue["RotationAnchorShift"],
anchorReference=OptionValue["RotationAnchorReference"],
numberOfDomains,preferredColours=OptionValue["Colours"],coloursToUse,colourTable,makeRotatedPolytope,polytopes,graphicList},

(* Check input *)
twoDimensionalQ=C===1;
If[(A*B*C)!=Length@domains,
Message[DomainPlot::InputMismatch];Abort[]];
If[AnyTrue[domains,NonNegative[#]\[Nand]IntegerQ[#]&],
Message[DomainPlot::InvalidDomainIndex];Abort[]];

(* Preparations *)
numberOfDomains=Max@domains;
coordinates=Flatten[Table[{i,j,k},
{i,0,A-1},{j,0,B-1},{k,0,C-1}],2];
If[twoDimensionalQ,coordinates=coordinates[[All,{1,2}]]];
coordinateDomainMap=Association@Thread[coordinates->domains];

coloursToUse=Join[preferredColours,
ColorData[96,#]&/@Range[numberOfDomains-Length@preferredColours]];
colourTable=Join[<|0->Transparent|>,Association@Thread[
Range@numberOfDomains->coloursToUse[[;;numberOfDomains]]]];

(* Fitting to given crystal system *)
M=InputCheck["GetCrystalFamilyMetric",
crystalFamily,If[twoDimensionalQ,"2D","3D"]];
coordinateDomainMap=Association@KeyValueMap[
M.#1->#2&,coordinateDomainMap];
If[graphicFunction=!=Automatic,
makePolytope=graphicFunction,
makePolytope[origin_]:=Parallelepiped[origin,Transpose@M]
];

(* Checking any rotations of domains *)
rotationQ=rotationMap=!=<||>;
If[rotationQ,
R=InputCheck["RotationTransformation",
{{A,B,C},domains},{anchorShift,anchorReference,rotationMap}];
makeRotatedPolytope[xyz_,d_]:=If[anchorReference==="Host"||(twoDimensionalQ&&anchorReference=!="Unit"),
GeometricTransformation[makePolytope[xyz],R[d]],
GeometricTransformation[makePolytope[xyz],R[d,xyz]]]
];

(* Preparing graphics *)
polytopes=If[rotationQ,
KeyValueMap[{colourTable@#2,makeRotatedPolytope[#1,#2]}&,coordinateDomainMap],
KeyValueMap[{colourTable@#2,makePolytope@#1}&,coordinateDomainMap]
];
graphicList={Opacity@OptionValue[Opacity],polytopes};

If[twoDimensionalQ,
Graphics[graphicList,Frame->False],
Graphics3D[graphicList,Boxed->False]
]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
EmbedStructure::InvalidGuestInput="Invalid guest unit input.";
EmbedStructure::InvalidTargetPositions="Invalid position input.";
EmbedStructure::InvalidProbabilities="The probabilities must be numbers between 0 and 1.";
EmbedStructure::InvalidTrimming="Invalid setting for \"TrimBoundary\".";
EmbedStructure::InvalidAlteration="Invalid input for \"Distortions\" or \"Rotations\".";
EmbedStructure::InvalidAlterationValues="Distortion/rotation amplitudes should be numeric and on the form \[Delta] or {\!\(\*SubscriptBox[\(\[Delta]\), \(min\)]\), \!\(\*SubscriptBox[\(\[Delta]\), \(max\)]\)}.";
EmbedStructure::InvalidDistortionType="\"DistortionType\" must be set to either \"Crystallographic\" or \"Cartesian\".";
EmbedStructure::VoidHost="Host structure cannot be 'Void'.";
EmbedStructure::InvalidOverlapRadius="\"OverlapRadius\" must be numeric.";
EmbedStructure::InvalidAnchorReference="Anchor reference type \[LeftGuillemet]`1`\[RightGuillemet] is invalid. Use either \"Host\" or \"Unit\".";

Options@EmbedStructure={
"DataFile"->FileNameJoin[{$MaXrdPath,"UserData","CrystalData.m"}],
"DistortionType"->"Cartesian",
"Distortions"->{0,0,0},
"MatchHostSize"->True,
"NewLabel"->"",
"OverlapPrecedence"->"",
"OverlapRadius"->1.0,
"RotationAnchorReference"->"Unit",
"RotationAnchorShift"->{0,0,0},
"RotationAxes"->IdentityMatrix[3],
"Rotations"->{0,0,0},
"ShowProgress"->False,
"TrimBoundary"->"None"
};

SyntaxInformation@EmbedStructure={
"ArgumentsPattern"->{_,_,_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
EmbedStructure[
guestUnitsInput_,
targetPositionsInput_List,
hostCrystal_String,
OptionsPattern[]
]:=Block[{
newStructureLabel=OptionValue["NewLabel"],
invAbort,conditionFilterQ=False,
crystalDataOriginal=$CrystalData,dataFile=OptionValue["DataFile"],
hostStructureSize,newSize,
probabilities,units,distributionList,i,
guestUnits,guestCopies,crystalLabels,nonVoidRange,
makeElementCrystal,
matchHostSizeQ=TrueQ@OptionValue["MatchHostSize"],
anchorShift=OptionValue["RotationAnchorShift"],
anchorReference=OptionValue["RotationAnchorReference"],
R,rotationAxes=OptionValue["RotationAxes"],
targetPositions=targetPositionsInput,embedLength,copyTranslations,hostCoordinates,mid,
latticeParameters,latticeParametersABC,hostM,hostMinverse,targetPositionsCartesian,
completed,M,T,p,P,CheckAndMakeRuleList,
distortions,rotations,performShift,performTwist,
conditions,list,shift,twist,
MakeAlteration,PrepareAlterationList,
conditionedDistortionsQ=False,conditionedRotationsQ=False,
coordinatesCrystal,coordinatesCartesian,
coordinatesCrystalEmbedded,coordinatesCrystalEmbeddedTranslated,
newCoordinates,newCoordinatesCartesian,
atomDataHost,atomDataGuest,joinedAtomData,boundary,hostCopy,newUnitCellAtomCount,
findOverlap,intervals,checks,atomData1,atomData2,xyz1,xyz2,nearestList,overlappingCoordinates,
overlapPrecedence=OptionValue["OverlapPrecedence"],
overlapRadius=OptionValue["OverlapRadius"]
},

(*--- Input checks ---*)
If[hostCrystal==="Void",Message[EmbedStructure::VoidHost];Abort[]];

boundary=OptionValue["TrimBoundary"];
If[!MemberQ[{"Box","None","OuterEdges"},boundary],
Message[EmbedStructure::InvalidTrimming];Abort[]];

If[!NumericQ@overlapRadius,
Message[EmbedStructure::InvalidOverlapRadius];Abort[]];

invAbort[]:=(Message[EmbedStructure::InvalidGuestInput];Abort[]);
Which[
(* A. 'guestUnits' as list of crystals or elements *)
ListQ@guestUnitsInput,
	If[guestUnitsInput==={},invAbort[]];
	Which[
	(* A.a. Regular crystal entries *)
	AllTrue[guestUnitsInput,StringQ]&&(Depth@guestUnitsInput===2),
	Null(* Check complete *),

	(* A.b. Conditional rules *)
	AllTrue[guestUnitsInput,Head[#]===Rule&]&&
	AllTrue[guestUnitsInput[[All,1]],Head[#]===Condition&],
	conditionFilterQ=True,

	True,invAbort[]
	],

(* B. 'guestUnits' as list paired with probabilities *)
Head@guestUnitsInput===Rule,
If[!MatchQ[Length/@guestUnitsInput,x_->x_/;x==x],invAbort[]];
If[!AllTrue[guestUnitsInput[[1]],0.0<=#<=1.0&],
Message[EmbedStructure::InvalidProbabilities];Abort[]],

(* Invalid input *)
True,invAbort[]
];

If[!MatchQ[Dimensions@targetPositionsInput,{_,3}],
Message[EmbedStructure::InvalidTargetPositions];Abort[]];

If[!MemberQ[{"Host","Unit"},anchorReference],
Message[EmbedStructure::InvalidAnchorReference,anchorReference];Abort[]];


(*--- Checking and preparing embeddings ---*)
crystalLabels=Cases[Flatten[
{hostCrystal,guestUnitsInput}],_String,3];
crystalLabels=DeleteCases[crystalLabels,"Void"];
makeElementCrystal[x_]:=
ImportCrystalData[
{x,x,"P1"},{1.,1.,1.,90.,90.,90.},{<|
"Element"->x,
"FractionalCoordinates"->{0.,0.,0.},
"DisplacementParameters"->0,
"Type"->"Uiso"|>},
"OverwriteWarning"->False];
makeElementCrystal/@Intersection[
crystalLabels,Keys@$PeriodicTable];
Scan[InputCheck["CrystalQ",#]&,crystalLabels];

overlapRadius=overlapRadius/GetLatticeParameters[
hostCrystal,"Units"->False][[{1,2,3}]];

hostCopy=$CrystalData[hostCrystal];
hostCoordinates=hostCopy[["AtomData",All,"FractionalCoordinates"]];
hostStructureSize=hostCopy["Notes"]["StructureSize"];
If[!ListQ@hostStructureSize,hostStructureSize={0,0,0}];

(*--- Preparing target positions ---*)
If[matchHostSizeQ&&hostStructureSize=!={0,0,0},
copyTranslations=Flatten[Table[
{i,j,k},{i,0,#1},{j,0,#2},{k,0,#3}]
&@@hostStructureSize,2];
targetPositions=Flatten[Outer[
Plus,copyTranslations,targetPositions,1],1];
targetPositions=DeleteCases[targetPositions,{x_,y_,z_}/;
Or@@MapThread[Greater,{{x,y,z},hostStructureSize}]];
(* If any negative coordinates, assume host is centred around origin *)
If[AnyTrue[Flatten@hostCoordinates,#<-1.&],
mid=\[LeftFloor]hostStructureSize/2.\[RightFloor];
targetPositions=#-mid&/@targetPositions
]];
embedLength=Length@targetPositions;

(* Preparing list to be used *)
guestUnits=Which[
conditionFilterQ,
targetPositions/.Append[guestUnitsInput,
{x_,y_,z_}/;True->"Void"],

Head@guestUnitsInput===Rule,
{probabilities,units}={Keys@#,Values@#}&@guestUnitsInput;
distributionList=Round[probabilities*embedLength];
Which[
Total@distributionList<embedLength,
While[Total@distributionList<embedLength,
i=RandomInteger[{1,Length@distributionList}];
distributionList[[i]]+=1],

Total@distributionList>embedLength,
While[Total@distributionList>embedLength,
i=RandomInteger[{1,Length@distributionList}];
distributionList[[i]]-=1]
];
RandomSample@Flatten@MapThread[
ConstantArray,{units,distributionList}],

True,
PadRight[#,embedLength,#]&@guestUnitsInput
];

guestCopies=$CrystalData/@guestUnits;
nonVoidRange=Complement[
Range@embedLength,
Flatten@Position[guestCopies,_Missing]];
If[nonVoidRange==={},Goto["End"]];

latticeParameters=GetLatticeParameters[
hostCrystal,"Units"->False];
latticeParametersABC=latticeParameters[[;;3]];
hostM=GetCrystalMetric[
hostCrystal,"ToCartesian"->True];
hostMinverse=Inverse@hostM;
targetPositionsCartesian=hostM.#&/@targetPositions;


(*--- Distortions and rotations -- Checks and preparations ---*)
MakeAlteration[c_]:=c;
MakeAlteration[{min_,max_}]:=Hold@RandomReal[{min,max}];
MakeAlteration[{x_,y_,z_}]:=MakeAlteration/@{x,y,z};
PrepareAlterationList[conditionsQ_,ruleList_List]:=If[conditionsQ,
conditions=Append[ruleList,{x_,y_,z_}/;True->{0.,0.,0.}];
list=Map[MakeAlteration,conditions,{2}];
#/.ReleaseHold@list&/@targetPositions,

ReleaseHold@ConstantArray[MakeAlteration/@ruleList,embedLength]
];

distortions=OptionValue["Distortions"];
performShift=distortions=!={0,0,0};

rotations=OptionValue["Rotations"];
performTwist=rotations=!={0,0,0};
rotations=rotations/.{
(c_Condition->r_List):>(N@c->N@r),
{r1_,r2_,r3_}:>N@{r1,r2,r3}
};

p=(NumericQ[#])&;P[x_]:=p[x];P[{x_,y_}]:=p[x]&&p[y];

CheckAndMakeRuleList[input_]:=Check[Which[
(input==={0,0,0})||SubsetQ[{Integer,Real,List},Head/@input],
	{False,CheckAndMakeRuleList[N@input,"Numeric"]},
(AllTrue[Input,Head[#]===Rule&]&&
AllTrue[input[[All,1]],Head[#]===Condition&])||(Head/@input==={Condition,List}),
	{True,CheckAndMakeRuleList[N@input,"Conditions"]},
True,
	Message[EmbedStructure::InvalidAlteration];
	Abort[]
],Abort[]];

	CheckAndMakeRuleList[input_,"Conditions"]:=(
	If[(P/@#)=!={True,True,True}&/@input[[All,2]],
	Message[EmbedStructure::InvalidAlterationValues];
	Abort[]];
	input/.{x_Condition->y_}:>{x->N@y}
	);

	CheckAndMakeRuleList[input_,"Numeric"]:=(
	If[(P/@input)=!={True,True,True},
	Message[EmbedStructure::InvalidAlterationValues];
	Abort[]];
	input
	);

{conditionedDistortionsQ,distortions}=CheckAndMakeRuleList@distortions;
distortions=PrepareAlterationList[conditionedDistortionsQ,distortions];
Which[
OptionValue["DistortionType"]==="Crystallographic",
Null,

OptionValue["DistortionType"]==="Cartesian",
distortions=hostMinverse.#&/@distortions,

True,
Message[EmbedStructure::InvalidDistortionType];Abort[]
];

{conditionedRotationsQ,rotations}=CheckAndMakeRuleList@rotations;
rotations=PrepareAlterationList[conditionedRotationsQ,rotations];

(*--- Actual tranformation -- Loop through each guest unit ---*)
R=If[performTwist,
InputCheck["RotationTransformation",{{0,0,0},{}},
{anchorShift,anchorReference,{},rotationAxes}]
];

completed=0;
If[TrueQ@OptionValue["ShowProgress"],
PrintTemporary[ProgressIndicator@Dynamic[completed/embedLength]]];

Do[
M=GetCrystalMetric[guestUnits[[i]],"ToCartesian"->True];
T=TranslationTransform[targetPositions[[i]]];

coordinatesCrystal=
guestCopies[[i,"AtomData",All,"FractionalCoordinates"]];

coordinatesCartesian=M.#&/@coordinatesCrystal;

coordinatesCrystalEmbedded=hostMinverse.#
&/@coordinatesCartesian;

coordinatesCrystalEmbeddedTranslated=T/@coordinatesCrystalEmbedded;
newCoordinates=coordinatesCrystalEmbeddedTranslated;

(* Optional: Perform shifts, twists or transformations *)
If[performShift||performTwist,

(* Optional: Rotations *)
If[performTwist,
twist=R[rotations[[i]],targetPositionsCartesian[[i]]];
newCoordinatesCartesian=hostM.#&/@newCoordinates;
newCoordinatesCartesian=twist@newCoordinatesCartesian;
newCoordinates=hostMinverse.#&/@newCoordinatesCartesian
];

(* Optional: Distortions *)
If[performShift,
shift=distortions[[i]];
newCoordinates=#+shift&/@newCoordinates]
];

guestCopies[[i,"AtomData",All,"FractionalCoordinates"]]=newCoordinates;
completed++,
{i,nonVoidRange}
];

(*--- Merge guest units with traget crystal ---*)
atomDataHost=$CrystalData[hostCrystal,"AtomData"];
atomDataGuest=Flatten@guestCopies[[nonVoidRange,"AtomData"]];

(* Optional: Remove overlapping target positions *)
joinedAtomData=If[!MemberQ[{"Host","Guest"},overlapPrecedence],
(* A. No measure taken against overlapping *)
Join[atomDataHost,atomDataGuest],

(* B. Remove superpositioned atoms *)
findOverlap[point_,region_List]:=(
If[region==={},Return@Nothing];
intervals=MapThread[Interval[{#1-#2,#1+#2}]&,{point,overlapRadius}];
checks=MapThread[IntervalMemberQ[#1,#2]&,{intervals,Transpose@region}];
checks=Position[And@@@Transpose@checks,True];
If[checks=!={},Part[region,Flatten@checks],Nothing]
);
{atomData1,atomData2}=If[overlapPrecedence==="Host",
#,Reverse@#]&@{atomDataHost,atomDataGuest};
{xyz1,xyz2}=Part[#,All,"FractionalCoordinates"]&/@{atomData1,atomData2};
nearestList=Nearest[xyz2,#,5]&/@xyz1;
overlappingCoordinates=Flatten[MapThread[
findOverlap,{xyz1,nearestList}],1];
atomData2=DeleteCases[atomData2,x_/;MemberQ[
overlappingCoordinates,x["FractionalCoordinates"]]];
Join[atomData1,atomData2]
];

(* Optional: Trim the outer boundary *)
If[boundary=!="None",
Which[
boundary==="Box",
joinedAtomData=DeleteCases[joinedAtomData,{x_,y_,z_}/;
Nand@@MapThread[0<=#1<#2&,{{x,y,z},hostStructureSize}],
{2}],

boundary==="OuterEdges",
joinedAtomData=DeleteCases[joinedAtomData,{x_,y_,z_}/;
Nand@@MapThread[0<=#1<#2&,{{x,y},hostStructureSize[[{1,2}]]}],
{2}]
];
joinedAtomData=Select[joinedAtomData,
KeyExistsQ[#,"FractionalCoordinates"]&]
];
joinedAtomData=MapAt[N@Chop[#,10.^-6]&,joinedAtomData,
{All,"FractionalCoordinates"}];
hostCopy["AtomData"]=joinedAtomData;

(* Overwrite host or create new crystal object *)
Label["End"];
If[newStructureLabel==="",newStructureLabel=hostCrystal];
$CrystalData=crystalDataOriginal;

newSize=\[LeftCeiling]targetPositions[[-1]]\[RightCeiling];
If[AnyTrue[newSize,#==0&],newSize+=1];
newUnitCellAtomCount=\[LeftCeiling]Length@joinedAtomData/Times@@newSize\[RightCeiling];
AppendTo[hostCopy["Notes"],"UnitCellAtomsCount"->newUnitCellAtomCount];
AppendTo[hostCopy["Notes"],"AsymmetricUnitAtomsCount"->Null];
AppendTo[hostCopy["Notes"],"StructureSize"->newSize];

InputCheck["Update$CrystalDataFile",dataFile,newStructureLabel,hostCopy];
newStructureLabel
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
SyntaxInformation@EquivalentIsotropicADP={
"ArgumentsPattern"->{_}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
EquivalentIsotropicADP[crystal_String]:=Block[{
calcEquiv,latticeParameters,abc\[Alpha]\[Beta]\[Gamma],a2b2c2,dispPars},

(* Input check *)
InputCheck["CrystalQ",crystal];

(* See:
  Fischer,R.X.& Tillmanns,E.(1988).Acta Cryst.C44,775-776.
https://doi.org/10.1107/S0108270187012745 *)
calcEquiv[
{a_,b_,c_,\[Alpha]_,\[Beta]_,\[Gamma]_},
{a2_,b2_,c2_},
{U11_,U22_,U33_,U12_,U13_,U23_}
]:=1/3*(
U11 (a*a2)^2+U22 (b*b2)^2+U33 (c*c2)^2+
2U12*a*b*a2*b2*Cos[\[Gamma]]+
2U13*a*c*a2*c2*Cos[\[Beta]]+
2U23*b*c*b2*c2*Cos[\[Alpha]]
);

calcEquiv[{__},{__},iso_]:=iso;

latticeParameters=GetLatticeParameters[
crystal,"Space"->"Both","Units"->False];
abc\[Alpha]\[Beta]\[Gamma]=latticeParameters[[1]]*{1,1,1,Degree,Degree,Degree};
a2b2c2=latticeParameters[[2,;;3]];
dispPars=$CrystalData[[crystal,"AtomData",All,"DisplacementParameters"]];

calcEquiv[abc\[Alpha]\[Beta]\[Gamma],a2b2c2,#]&/@dispPars
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
ExpandCrystal::InvalidSize="The structure size must be a list of three natural numbers.";
ExpandCrystal::DuplicateLabel="The new label must be different from the input.";

Options@ExpandCrystal={
"DataFile"->FileNameJoin[{$MaXrdPath,"UserData","CrystalData.m"}],
"ExpandIntoNegative"->False,
"FirstTransformTo"->False,
"IncludeBoundary"->True,
"NewLabel"->"",
"StoreTemporarily"->False
};

SyntaxInformation@ExpandCrystal={
"ArgumentsPattern"->{_,_.,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
ExpandCrystal[crystal_String,structureSize_List:{1,1,1},
OptionsPattern[]]:=Block[{
crystalDataOriginal=$CrystalData,dataFile=OptionValue["DataFile"],
newLabel=OptionValue["NewLabel"],
changeCell=OptionValue["FirstTransformTo"],
storeTempQ=TrueQ@OptionValue["StoreTemporarily"],
crystalData,crystalCopy,atomData,coordinates,spaceGroup,generated,
copyTranslations,mid,atomDataMapUnitCell,
cutoffFunction,atomDataMapExpanded,lengths,
newAtomData
},

(* Input check and data acquisition *)
If[
AllTrue[structureSize,Positive[#]&&IntegerQ[#]&]\[Nand]
Length[structureSize]===3,
Message[ExpandCrystal::InvalidSize];Abort[]];

InputCheck["CrystalQ",crystal];
crystalData=crystalCopy=$CrystalData[crystal];

If[crystal===newLabel,
Message[ExpandCrystal::DuplicateLabel,newLabel];Abort[]];
If[newLabel==="",newLabel=crystal<>"_"<>StringJoin@Riffle[
ToString/@structureSize,"x"]];

(* Optional: Transform cell beforehand *)
If[changeCell=!=False,
InputCheck["InterpretSpaceGroup",changeCell];
AssociateTo[$CrystalData,newLabel->crystalCopy];
UnitCellTransformation[newLabel,changeCell];
crystalData=crystalCopy=$CrystalData[newLabel]
];

(* Expand asymmetric unit to unit cell *)
atomData=crystalData["AtomData"];
coordinates=atomData[[All,"FractionalCoordinates"]];
spaceGroup=crystalData["SpaceGroup"];
generated=N[SymmetryEquivalentPositions[
spaceGroup,#]&/@coordinates];

(* Generate full content of the unit cell *)
atomDataMapUnitCell=Association@Thread[
Range@Length@atomData->generated];

(* Copy by translation *)
copyTranslations=Flatten[
Table[{i,j,k},{i,0,#1},{j,0,#2},{k,0,#3}]&@@N@structureSize,
2];

atomDataMapExpanded=Flatten[
Outer[Plus,copyTranslations,#,1],
1]&/@atomDataMapUnitCell;

(* Optional: Complete the outer boundary *)
cutoffFunction=If[TrueQ@OptionValue["IncludeBoundary"],
Greater,GreaterEqual];

(* Delete atoms whose coordinates are outside *)
atomDataMapExpanded=DeleteCases[atomDataMapExpanded,{x_,y_,z_}/;
Or@@MapThread[cutoffFunction,{{x,y,z},structureSize}],
{2}];

(* Optional: Center translations around origin *)
If[TrueQ@OptionValue["ExpandIntoNegative"],
mid=\[LeftFloor]structureSize/2.\[RightFloor];
atomDataMapExpanded=Map[#-mid&,atomDataMapExpanded,{2}];
];

(* Create new atom data structure *)
lengths=Values[Length/@atomDataMapExpanded];
newAtomData=Table[
ConstantArray[atomData[[i]],lengths[[i]]],
{i,Length@atomData}];

newAtomData[[All,All,"FractionalCoordinates"]]=
Values@atomDataMapExpanded;
newAtomData=Flatten[newAtomData,1];

(* Create new crystal entry *)
crystalCopy["AtomData"]=newAtomData;
AssociateTo[crystalCopy,"Notes"-><|
"StructureSize"->structureSize,
"UnitCellAtomsCount"->Total[Length/@generated],
"AsymmetricUnitAtomsCount"->Length@atomDataMapUnitCell
|>];

(* If temporary storage was used, reset pointer to original $CrystalData *)
If[storeTempQ,
MaXrd`Private`$TempCrystalData=<|newLabel->crystalCopy|>;
$CrystalData=crystalDataOriginal;
InputCheck["Update$CrystalDataAutoCompletion"],

InputCheck["Update$CrystalDataFile",dataFile,newLabel,crystalCopy]
];

newLabel
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
ExportCrystalData::InvalidProgramOrFormat="\[LeftGuillemet]`1`\[RightGuillemet] is not a valid program/format.";
ExportCrystalData::ParameterError="Invalid input parameters.";
ExportCrystalData::DirectoryExpected="An existing output directory was expected.";
ExportCrystalData::InvalidSubtractionMode="\[LeftGuillemet]`1`\[RightGuillemet] is not a valid scattering subtraction mode.";

Options@ExportCrystalData={
"Detailed"->False,
"IncludeStructureSizeInfo"->True
};

SyntaxInformation@ExportCrystalData={
"ArgumentsPattern"->{_,_,_,_.,_.,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
ExportCrystalData["DISCUS",crystal_String,outputFile_String,OptionsPattern[]]:=Block[{
crystalData,atomData,size,unitCellAtomCount,latticeParameters,
formatter,appendComma,simpleQ,
preambleTitle,preambleSpaceGroup,preambleCell,preambleCount,preambleAtomsHeader,preamble,
elements,coordinates,dispPars,items,spacesToUse,additional,makeSpace,atoms
},

(* Loading necessary data *)
{crystalData,atomData,size,latticeParameters}=ECD$LoadNecessaries@crystal;
unitCellAtomCount=Lookup[crystalData["Notes"],
"UnitCellAtomsCount",Round[Length[atomData]/Times@@size]];

(* Auxiliary *)
formatter[x_]:=ToString@NumberForm[
N@Chop[x,Power[10.,-5]],{9,6}];
appendComma[x_]:=Map[StringInsert[#,",",-1]&,x,{1}];
simpleQ=!TrueQ@OptionValue["Detailed"];

latticeParameters=If[simpleQ,
ToString/@N@latticeParameters,
Map[formatter,latticeParameters,{1}]];

(* Creating the preamble *)
preambleTitle={"title  ",crystal};
preambleSpaceGroup={"spcgr  ",StringDelete[
$GroupSymbolRedirect[crystalData["SpaceGroup"]]
["Name","HermannMauguinShort"]," "]};
preambleCell={"cell   ",StringJoin@Riffle[latticeParameters,"  "]};
preambleCount={"ncell  ",StringJoin@Riffle[
ToString/@Join[size,{unitCellAtomCount}],
",  "]};
preambleAtomsHeader=If[simpleQ,
{"atoms\n"},
{"atoms      x,              y,              z,             Biso,    Property,  MoleNo,  MoleAt,   Occ\n"}];
preamble=StringJoin@Riffle[{
preambleTitle,preambleSpaceGroup,preambleCell,
If[TrueQ@OptionValue["IncludeStructureSizeInfo"],
preambleCount,Nothing],
preambleAtomsHeader},
{"\n"}
];

(* Extracting atom data *)
elements=ToUpperCase@Lookup[atomData,"Element"];
elements=StringDelete[elements,{"+","-",DigitCharacter}];
coordinates=N@Lookup[atomData,"FractionalCoordinates"];
coordinates=Map[formatter,coordinates,{2}];
If[!simpleQ,coordinates=appendComma@coordinates];
dispPars=EquivalentIsotropicADP[crystal]/._Missing->0.;
dispPars=formatter/@dispPars;
If[!simpleQ,dispPars=appendComma@dispPars];
items={elements,Sequence@@Transpose@coordinates,dispPars};

(* Determining spacing information *)
spacesToUse=Transpose@ConstantArray[
{11,16,16,15,8},Length@atomData];
additional="1,       0,       0,   1.000000";
Do[
spacesToUse[[i]]=spacesToUse[[i]]-(StringLength/@items[[i]]),
{i,Length@items-1}];
spacesToUse=Transpose@spacesToUse;
makeSpace[i_,j_]:=ConstantArray[" ",spacesToUse[[i,j]]];

(* Writing atom data *)
atoms=Reap[Do[Sow[{
elements[[i]],makeSpace[i,1],
coordinates[[i,1]],makeSpace[i,2],
coordinates[[i,2]],makeSpace[i,3],
coordinates[[i,3]],makeSpace[i,4],
dispPars[[i]],makeSpace[i,5],
If[simpleQ,"",additional],
"\n"
}],{i,Length@atomData}]][[2,1]];

(* Prepare output and export *)
Export[outputFile,StringJoin[preamble,atoms],"String"]
]


(* ::Input::Initialization:: *)
ExportCrystalData["DIFFUSE",crystal_String,outputDir_String,hklPlane_,indexLimit_,subtractionMode_String,OptionsPattern[]]:=Block[{
reorganise,source,scatteringFactorTemplate,
crystalData,atomData,allElements,size,latticeParameters,directionCosines,M,
partA,X,Y,Z,unitCells,MakeSitesTable,table,sitesList,
partB,maxPossibleSites,newAtomData,blockLengths,
headerComments,headerData,padWidth,header,scatteringData
},

(* Checks *)
InputCheck["CrystalQ",crystal];
If[!DirectoryQ@outputDir,
Message[ExportCrystalData::DirectoryExpected];Abort[]];
If[!MemberQ[{"N","Y","E","e"},subtractionMode],
Message[ExportCrystalData::InvalidSubtractionMode,subtractionMode];Abort[]];

(* Auxiliary *)
reorganise[data_Association]:=With[{xyz=data["FractionalCoordinates"]},
{Min/@Transpose[{IntegerPart@xyz+1,size}],
StringJoin[
"   ",
StringPadRight[ToString@DecimalForm[#,{11,5}]&/@FractionalPart@xyz,11],
StringDelete[data["Element"],{DigitCharacter,"+","-"}]
]}
];

MakeSitesTable[stop_Integer]:=(
table=Table[{i,If[i>stop,0,i]},{i,maxPossibleSites}];
table=Map[ToString,table,{2}];
table=StringPadLeft[#,6]&/@table;
StringJoin/@table
);

source=Import@FileNameJoin[
{$MaXrdPath,"Core","Data","AtomicScatteringFactor","InternationalTablesC(3rd).m"}];
scatteringFactorTemplate[element_]:=StringTemplate["'`X`'
`a1`,`b1`,`a2`,`b2`
`a3`,`b3`,`a4`,`b4`,`c`
0.0,0.0"]@Join[source[element],<|"X"->element|>];

(* Loading necessary data *)
{crystalData,atomData,size,latticeParameters}=ECD$LoadNecessaries@crystal;
allElements=DeleteDuplicates@atomData[[All,"Element"]];
allElements=DeleteDuplicates[StringDelete[#,
{DigitCharacter,"+","-"}]&/@allElements];
directionCosines=ToString@DecimalForm[
Cos[#*Degree],{6,5}]&/@N@latticeParameters[[4;;6]];
M=GetCrystalMetric[crystal,"ToCartesian"->True];

(* File 1: Crystal data in 'DIFFUSE' format *)
(* Part B *)
newAtomData=reorganise/@atomData;
newAtomData=DeleteDuplicates/@GatherBy[newAtomData,#[[1]]&];
blockLengths=Length/@newAtomData;
maxPossibleSites=Max@blockLengths;
partB=Flatten@Riffle[
StringTemplate[" `1` `2` `3`"]@@@newAtomData[[All,1,1]],
newAtomData[[All,All,2]]
];

(* Part A *)
sitesList=MakeSitesTable/@blockLengths;
{X,Y,Z}=size;
unitCells=Flatten[Table[
StringTemplate[" `1` `2` `3`"][x,y,z],
{x,1,X},{y,1,Y},{z,1,Z}],2];
partA=Flatten@Transpose[{unitCells,sitesList}];
PrependTo[partA,StringTemplate[" `X` `Y` `Z` `N1` `N1`"]@<|
"X"->X,"Y"->Y,"Z"->Z,"N1"->maxPossibleSites|>];

(* File 2: Input/summary file *)
headerComments={
"Header or structure label",
"Random number seeds",
"Cell coord's; angles are cos(ang)",
"Size of crystal simulation (unit cells)",
"Periodic Boundary?",
"Origin of volume",
"u-axis (bottom right) and image x-dimension",
"v-axis (top left) and image y-dimension",
"w-axis (top right) and image z-dimension",
"sin(theta)/lambda maximum",
"Lot size",
"Number of lots",
"Number of atom sites per cell",
"Number of different atom types (list of sc.coef's below)",
"Subtract average lattice? ('N', 'e', 'E' or 'Y')"
};
headerData={
crystal,
"12645, 9676",
StringTemplate["`1` `2` `3`  `4` `5` `6`"]@@Join[latticeParameters[[1;;3]],directionCosines],
StringTemplate["`1` `2` `3`"]@@size,
"Y",
Sequence@@InputCheck["GetReciprocalImageOrientation",
crystal,hklPlane,indexLimit,True],
"3.0",
StringTemplate["`1` `2` `3`"]@@size,
"1",
ToString@maxPossibleSites,
ToString@Length@allElements,
subtractionMode
};
padWidth=Max[StringLength/@headerData]+3;
headerData=StringPadRight[headerData,padWidth];
header=Transpose@{headerData,"! "<>#&/@headerComments};
scatteringData=scatteringFactorTemplate/@allElements;

(* Prepare output and export *)
{
Export[
FileNameJoin[{outputDir,"diffuse_input1_crystal.txt"}],
Join[partA,partB,{"\n"}],
"Table"],
Export[
FileNameJoin[{outputDir,"diffuse_input2_setup.txt"}],
Join[header,scatteringData,{"\n"}],
"Table"]
}
]


(* ::Input::Initialization:: *)
ECD$LoadNecessaries[crystal_String]:=Block[
{crystalData,atomData,crystalNotes,size,latticeParameters},
InputCheck["CrystalQ",crystal];
crystalData=$CrystalData[crystal];
atomData=N@crystalData["AtomData"];
crystalNotes=Lookup[crystalData,"Notes",<||>];
If[!AssociationQ@crystalNotes,crystalNotes=<||>];
size=Lookup[crystalNotes,"StructureSize",
Round/@Max/@Transpose@atomData[[All,"FractionalCoordinates"]]
]/.{0->1};
latticeParameters=GetLatticeParameters[crystal,"Units"->False];
{crystalData,atomData,size,latticeParameters}
]


(* ::Input::Initialization:: *)
ExportCrystalData[invalidProgram_String,rest___]:=(
If[!MemberQ[{"DISCUS","DIFFUSE"},invalidProgram],
Message[ExportCrystalData::InvalidProgramOrFormat,invalidProgram],
Message[ExportCrystalData::ParameterError]];
Abort[]
)


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
Options@ExtinctionLength={
"Polarisation"->"\[Pi]",
"Units"->True
};

SyntaxInformation@ExtinctionLength={
"ArgumentsPattern"->{_,_,_.,_.,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
ExtinctionLength[
crystal_,
lambda:_?(NumericQ[#]||QuantityQ[#]&):-1,
hklInput_List,
\[Gamma]o:_?NumericQ:1,
\[Gamma]h:_?NumericQ:1,
OptionsPattern[]]:=Block[{
sg,hkl,L,
\[Lambda],V,\[Theta],R,C,FhFhbar,g,\[CapitalLambda]o,
temp},

(* Check input *)
InputCheck["CrystalQ",crystal];
sg=$CrystalData[crystal,"SpaceGroup"];
hkl=InputCheck[hklInput,"Integer","WrapSingle"];
L=Length@hkl;

(* Wavelength *)
\[Lambda]=InputCheck["ProcessWavelength",crystal,lambda];

(* Unit cell volume *)
V=Sqrt@Det@GetCrystalMetric[crystal];

(* Classical electron radius *)
R=If[OptionValue["Units"],
QuantityMagnitude@UnitConvert[
Quantity["ClassicalElectronRadius"],"Angstroms"],
2.81794032*Power[10,-5](* \[ARing]ngstr\[ODoubleDot]ms *)
];

(* Bragg angle and polarisation *)
\[Theta]=BraggAngle[crystal,\[Lambda],hkl,"Units"->False]*Degree;
C=Flatten[{InputCheck["Polarisation",
OptionValue["Polarisation"],2\[Theta]]}];

(* Structure factors *)
temp=StructureFactor[crystal,#,\[Lambda],"Units"->False]&/@{hkl,-hkl};
temp=Cases[temp,{F_?NumericQ,\[Phi]_?NumericQ}:>F,3];
FhFhbar=Times@@ArrayReshape[temp,{2,L}];

	(* Message about extinction *)
	Do[If[FhFhbar[[i]]==0,
	Message[StructureFactor::extinct,hkl[[i]],sg]
	],{i,L}];

(* Extinction (Pendell\[ODoubleDot]sung) distance *)
If[\[Gamma]o===\[Gamma]h===1,g=1,g=\[Pi]*Sqrt[\[Gamma]o*Abs[\[Gamma]h]]];(* geometrical factor *)
\[CapitalLambda]o=Reap[Do[Quiet@Sow[
(V *g)/(R*\[Lambda]*C[[i]]*Sqrt@FhFhbar[[i]])*(1 (* \[Mu]m *)/(10^4) (* \[CapitalARing] *))],
{i,L}]][[2,1]];
\[CapitalLambda]o=\[CapitalLambda]o/.ComplexInfinity->Undefined;

(* Option: Units *)
If[OptionValue["Units"],
\[CapitalLambda]o=\[CapitalLambda]o/.x_?NumericQ:>Quantity[x,"Micrometers"]];

If[L===1,First@\[CapitalLambda]o,\[CapitalLambda]o]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
GetAtomicScatteringFactors::source="Invalid data source option.";
GetAtomicScatteringFactors::missing="\[LeftGuillemet]`1`\[RightGuillemet] is missing in data source for the `2`.";
GetAtomicScatteringFactors::slRequired="Crystal name or Sin[\[Theta]]/\[Lambda] values required.";
GetAtomicScatteringFactors::slRange="The value `1` \!\(\*SuperscriptBox[\(\[CapitalARing]\), \(-1\)]\) is out of range for the f0 source `2`.";

Options@GetAtomicScatteringFactors={
"DispersionCorrections"->True,
"f0Source"->"WaasmaierKirfel",
"f1f2Source"->"CromerLiberman",
"IgnoreIonCharge"->True,
"SeparateCorrections"->False
};

SyntaxInformation@GetAtomicScatteringFactors={
"ArgumentsPattern"->{_,_,_.,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
GetAtomicScatteringFactors[
crystal_String,
hklInput_List,
input\[Lambda]:_?(NumericQ[#]||QuantityQ[#]&):-1,
OptionsPattern[]]:=Block[{
\[Lambda]=input\[Lambda],hkl,elements,sl,options},

(*---* Basic *---*)
(* Crystal and wavelength/energy *)
\[Lambda]=InputCheck["ProcessWavelength",crystal,\[Lambda],False];

(* Reflection(s) *)
hkl=InputCheck[hklInput,"WrapSingle","Integer"];

(* Processing elements *)
elements=GetElements[crystal,"Tally"->False];

	(* Optional: Remove charge of ions *)
	If[OptionValue["IgnoreIonCharge"],
	elements=DeleteDuplicates@
	StringDelete[elements,
	{DigitCharacter,"+","-"}]];

(* Sin[\[Theta]]/\[Lambda] *)
H=Chop@N@Inverse@GetCrystalMetric[crystal];
sl=Sqrt[#.H.#]/2.&/@hkl;

(*---* Relaying data to main function *---*)
options=#->OptionValue[#]&/@
Keys@Options@GetAtomicScatteringFactors;

GetAtomicScatteringFactors[elements,sl,\[Lambda],Sequence@@options]
]


(* ::Input::Initialization:: *)
GetAtomicScatteringFactors[
inputElements_List,
inputSL_List,
input\[Lambda]:_?(NumericQ[#]||QuantityQ[#]&):-1,
OptionsPattern[]]:=Block[{
\[Lambda]=input\[Lambda],elements=inputElements,sl=inputSL,
f0Source,f1f2Source,$f0Local,$f1f2Local,
upperLimit,ignore,
addCorrectionsQ=TrueQ@OptionValue["DispersionCorrections"],
separateQ=TrueQ@OptionValue["SeparateCorrections"],
ipfQ,
coefficients,akeys,bkeys,a,b,c,
f0,corrections,output,
temp},
(*---* Checking input *---*)
(* Elements *)
elements=InputCheck["InterpretElement",Flatten[{elements}]];

(* Sin[\[Theta]]/\[Lambda] *)
sl=Flatten[{sl}];
If[!AllTrue[sl,TrueQ@Not@Negative[#]&],
Message[GetAtomicScatteringFactors::slRequired];Abort[]];

(* Wavelength *)
If[\[Lambda]!=-1,\[Lambda]=QuantityMagnitude@
InputCheck["GetEnergyWavelength",\[Lambda]]];	

(* Data sources *)
f0Source=OptionValue["f0Source"];
f1f2Source=OptionValue["f1f2Source"];
ignore={"H","He"};

	(* Validating sources *)
	If[!MemberQ[FileBaseName/@FileNames["*.m",
FileNameJoin[{
$MaXrdPath,"Core","Data",
"AtomicScatteringFactor"}]],
f0Source]||
	!MemberQ[FileBaseName/@FileNames["*.m",
FileNameJoin[{
$MaXrdPath,"Core","Data",
"AtomicScatteringFactor","AnomalousCorrections"}]],
f1f2Source],
	Message[GetAtomicScatteringFactors::source];Abort[]];

(*---* Useful variables *---*)
(* Check specific range limits *)
upperLimit=Which[
f0Source==="CromerMann",1.5,
f0Source==="InternationalTablesC(3rd)",2.0,
f0Source==="WaasmaierKirfel",6.0,
True,2.5
];

Do[
temp=sl[[i]];
If[temp>upperLimit,
Message[GetAtomicScatteringFactors::slRange,
ToString@temp,f0Source];
Abort[]],
{i,Length@sl}];

(* Loading source for calculating f0 *)
	(* Setup accumulative variable for the session *)
	If[!AssociationQ@$f0,$f0=<||>];
	If[
	(* a. Check if same source is imported already *)
	KeyExistsQ[$f0,f0Source],
	$f0Local=$f0[f0Source],

	(* b. Import specified data *)
	$f0Local=Import@FileNameJoin[{
	$MaXrdPath,"Core","Data",
	"AtomicScatteringFactor",f0Source<>".m"}];
		(* Update the accumulative variable *)
		AppendTo[$f0,f0Source->$f0Local]
	];

(* Check if atom types are found in $f0 source *)
temp=Complement[elements,Keys@$f0Local];
If[temp=!={},
Message[GetAtomicScatteringFactors::missing,
First@temp,"non-dispersive part (\!\(\*FormBox[SubscriptBox[\(f\), \(0\)],
TraditionalForm]\))"];
Abort[]];

(*---* Calculating form factor (f0) from tabulated data *---*)
ipfQ=Head@First@$f0Local===InterpolatingFunction;
Which[
(* a. Interpolation data *)
(* sin(\[Theta])/\[Lambda] *)
ipfQ,
	(* Non-dispersive part of form factor *)
	f0=Table[{X,$f0Local[X][s]},
	{s,sl},{X,elements}],

(* b. Coefficients *)
True,
	(* Stored with alternating 'a' and 'b'
	and 'c' last *)		
	coefficients=$f0Local[[elements]];

	{akeys,bkeys}=Flatten@StringCases[
Keys@First@coefficients,
#~~DigitCharacter..]&/@{"a","b"};

	{a,b,c}={
	Values/@coefficients[[All,akeys]],
	Values/@coefficients[[All,bkeys]],
	Values@coefficients[[All,"c"]]};

	(* Non-dispersive part of form factor *)
	f0=Table[
	{Keys[coefficients][[i]],
	Total[a[[i]]*Exp[-b[[i]]*(sl[[j]])^2]]+c[[i]]},
	{j,Length@sl},
	{i,Length@elements}]
];

(* Check: Correct normalisation by electrons *)
If[
Abs[If[ipfQ,
$f0Local["C"][0],
Total@Values@$f0Local[["C",akeys]]
+$f0Local["C","c"]]-6(* Carbon: Z=6 *)]>0.5,

f0=f0/.{X_String,f_?NumericQ}:>{X,f*$PeriodicTable[X,"AtomicNumber"]}
];

(*---* Dispersion corrections (f1 + f2) *---*)
	(* Optional: No dispersion corrections *)
	If[\[Lambda]===-1,addCorrectionsQ=False];
	If[!addCorrectionsQ,
	(* Prepare output *)
	output=Map[Association,MapThread[Rule,Transpose[#]]&/@f0];
	Goto["End"]];

	(* Loading source for calculating f' + f'' *)
		(* Setup global variable for the session *)
		If[!AssociationQ@$f1f2,$f1f2=<||>];
		If[
		(* a. Check if same source is imported already *)
		KeyExistsQ[$f1f2,f1f2Source],
		$f1f2Local=$f1f2[f1f2Source],

		(* b. Import specified data *)
		$f1f2Local=Import@FileNameJoin[{
	$MaXrdPath,"Core","Data",
	"AtomicScatteringFactor","AnomalousCorrections",
	f1f2Source<>".m"}];
			(* Update the session's global variable *)
			AppendTo[$f1f2,f1f2Source->$f1f2Local]
		];

	(* Check if atom types are found in $f1f2 source *)
	temp=Complement[elements,Keys@$f1f2Local];
	temp=temp/.Thread[ignore->Nothing];(* ignored elements *)
	If[temp=!={},
	Message[GetAtomicScatteringFactors::missing,
	First@temp,"dispersion corrections (\!\(\*FormBox[\(\*SuperscriptBox[\"f\", \"\[Prime]\",\nMultilineFunction->None] + \*SuperscriptBox[\"f\", \"\[Prime]\[Prime]\",\nMultilineFunction->None]\),
TraditionalForm]\))"];
	Abort[]];

	(* Procedure *)
	Do[
	If[MemberQ[ignore,elements[[j]]],Continue[]];
	corrections=$f1f2Local[elements[[j]]][\[Lambda]];
	AppendTo[f0[[i,j]],corrections],
		{i,Length@sl},
		{j,Length@elements}];

(*---* Preparing and returning association *---*)
If[separateQ,
(* a. Return f0 and f1f2 separated *)
f0=f0/.{
(* Elements with corrections *)
{X_String,f0_?NumericQ,f1f2_?NumericQ}:>
(X-><|"f0"->f0,"f1f2"->f1f2|>),
(* Without corrections (ignored elements) *)
{X_String,f0_?NumericQ}:>(X-><|"f0"->f0,"f1f2"->0.|>)
},

(* b. Return f0 and f1f2 combined *)
f0=f0/.{
(* Elements with corrections *)
{X_String,f0_?NumericQ,f1f2_?NumericQ}:>
(X->f0+f1f2),
(* Without corections (ignored elements) *)
{X_String,f0_?NumericQ}:>(X->f0)
}
];
output=Association/@f0;

Label["End"];
If[Length@output==1,First@output,output]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
GetCrystalMetric::InvalidInput="Invalid input.";
GetCrystalMetric::InvalidSpace="\"Space\" must either be \"Direct\" or \"Reciprocal\".";

Options@GetCrystalMetric={
"Space"->"Direct",
"ToCartesian"->False
};

SyntaxInformation@GetCrystalMetric={
"ArgumentsPattern"->{_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
GetCrystalMetric[userInput_,OptionsPattern[]]:=Block[{
input=userInput,
a,b,c,\[Alpha],\[Beta],\[Gamma],
c1,c2,c3,M,
space=OptionValue["Space"],
MakeMetric,ExtractParametersFromMatrix
},

(* Auxiliary *)
MakeMetric[{a_,b_,c_,alpha_,beta_,gamma_}]:=(
{\[Alpha],\[Beta],\[Gamma]}={alpha,beta,gamma};
If[AnyTrue[{\[Alpha],\[Beta],\[Gamma]},#>2\[Pi]&],{\[Alpha],\[Beta],\[Gamma]}*=Degree];
N@Chop[{
{a^2,a*b*Cos[\[Gamma]],a*c*Cos[\[Beta]]},
{a*b*Cos[\[Gamma]],b^2,b*c*Cos[\[Alpha]]},
{a*c*Cos[\[Beta]],b*c*Cos[\[Alpha]],c^2}}]
);

ExtractParametersFromMatrix[matrix_]:=(
{a,b,c}=Sqrt@Diagonal@matrix;
\[Alpha]=ArcCos[matrix[[2,3]]/(b*c)]/Degree;
\[Beta]=ArcCos[matrix[[1,3]]/(a*c)]/Degree;
\[Gamma]=ArcCos[matrix[[1,2]]/(a*b)]/Degree;
{a,b,c,\[Alpha],\[Beta],\[Gamma]}
);

(* Input check *)
If[!MemberQ[{"Direct","Reciprocal"},space],
Message[GetCrystalMetric::InvalidSpace];Abort[]
];

If[
StringQ@userInput,
(* A. Crystal label *)
InputCheck["CrystalQ",userInput];
{a,b,c,\[Alpha],\[Beta],\[Gamma]}=GetLatticeParameters[userInput,
"Space"->space,"Units"->False],

(* B. Lattice parameters directly *)
If[AssociationQ@input,
input=Lookup[userInput,{"a","b","c","\[Alpha]","\[Beta]","\[Gamma]"}]];

If[
!AllTrue[input,QuantityQ[#]||NumericQ[#]&]||
Length@Flatten@input!=6,
Message[GetCrystalMetric::InvalidInput];Abort[]];

{a,b,c,\[Alpha],\[Beta],\[Gamma]}=N@input;
If[AllTrue[input,QuantityQ],
{a,b,c}=QuantityMagnitude@UnitConvert[{a,b,c},"Angstroms"];
{\[Alpha],\[Beta],\[Gamma]}=QuantityMagnitude@UnitConvert[{\[Alpha],\[Beta],\[Gamma]},"Degrees"]
];

(* Optional: Use metric for reciprocal space *)
If[OptionValue["Space"]==="Reciprocal",
M=MakeMetric[{a,b,c,\[Alpha],\[Beta],\[Gamma]}];
{a,b,c,\[Alpha],\[Beta],\[Gamma]}=ExtractParametersFromMatrix@Inverse@M]
];
If[AnyTrue[{\[Alpha],\[Beta],\[Gamma]},#>2.\[Pi]&],{\[Alpha],\[Beta],\[Gamma]}*=Degree];

(* Metric tensor *)
M=If[TrueQ@OptionValue["ToCartesian"],
(* Optional: Return matrix that converts to Cartesian coordinates *)
N@Chop[{
{a,b*Cos[\[Gamma]],c1},
{0.,b*Sin[\[Gamma]],c2},
{0.,0.,c3}
}//.{
c1->c*Cos[\[Beta]],
c2->c*(Cos[\[Alpha]]-Cos[\[Gamma]]*Cos[\[Beta]])/Sin[\[Gamma]],
c3->Sqrt[c^2-c1^2-c2^2]
}],

MakeMetric[{a,b,c,\[Alpha],\[Beta],\[Gamma]}]
];

M
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
GetElements::InvalidFormula="Invalid chemical formula.";
GetElements::InvalidElements="Invalid elements detected: `1`.";

Options@GetElements={
"IgnoreIonCharge"->True,
"Tally"->False
};

SetAttributes[GetElements,Listable];

SyntaxInformation@GetElements={
"ArgumentsPattern"->{_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
GetElements[input_String,OptionsPattern[]]:=Block[{formula=input,patternX,groupX,tallyQ=TrueQ@OptionValue["Tally"],elements,temp},
(*---* Check input string *---*)
(*--* A. Crystal name input *--*)
If[MemberQ[Keys@$CrystalData,input],
(* a. Use chemical formula *)
formula=Lookup[$CrystalData@input,"ChemicalFormula",

(* b. Return elements contained in 'AtomData' *)
elements=$CrystalData[[input,"AtomData",All,"Element"]];
If[tallyQ,
(* Return tally *)
Return@Tally@elements,
(* Return elements only *)
Return@DeleteDuplicates@elements
]]];

(*--* B. Chemical formula string *--*)
If[StringContainsQ[ToString@FullForm@formula,"!"],
	(* a. Formatted string *)
		(* Considering the full form *)
		elements=ToString@FullForm@formula;

		(* Marking subscripts with '$' 
		elements=StringReplace[elements,"\\), \\("\[Rule]"$"];*)

		(* Cleaning the full form *)
		elements=StringDelete[elements,
		{"\\!\\(\\*SubscriptBox[\\(",
		"[\\(","\\)]\\)","\\), \\("}],

	(* b. Plain string *)
		elements=formula];

(*---* Useful local variables *---*)
	patternX={_?UpperCaseQ~~_?LowerCaseQ,_?UpperCaseQ};
	groupX[S_]:=StringCases[S,{
	x:patternX~~{n1:DigitCharacter..~~"."~~n2:DigitCharacter..}
	:>{x,ToExpression[n1<>"."<>n2]},
	x:patternX~~n:DigitCharacter..~~pm:{"+","-"}:>{x,n~~pm},
	x:patternX~~pm:{"+","-"}~~n:DigitCharacter...:>{x,n~~pm},
	x:patternX~~n:DigitCharacter...:>{x,ToExpression@n}
	}]/.{""->"1",Null->1};

(*---* Extracting symbols and numbers *---*)
	(* Distribute parenthesis subscripts *)
	elements=StringReplace[elements,
	"("~~p__~~")"~~s:DigitCharacter..
		:>StringJoin[ToString/@Flatten@MapAt[
#*ToExpression[s]&,groupX@p,{All,2}]]
	];

	(* Group elements and corresponding subscripts *)
	elements=groupX@elements;

		(* Check *)
		If[elements==={},
		Message[GetElements::InvalidFormula];Abort[]];
		
		temp=elements[[All,1]];
		temp=InputCheck["InterpretElement",temp];
		temp=DeleteDuplicates@temp;
		temp=Complement[temp,Keys@$PeriodicTable];
		If[temp=!={},
		Message[GetElements::InvalidElements,ToString@temp];
		Abort[]];

	(* Merge equal elements *)
	If[!DuplicateFreeQ[First/@elements],
	elements=MapAt[ToExpression,elements,{All,2}];
	elements=GatherBy[elements,First];
	elements=elements/.x_/;Depth[x]===3:>
	{x[[1,1]],Total@x[[All,2]]}];

	(* Checking for ions *)
	elements=elements/.x_List/;
Depth[x]==2&&StringContainsQ[
ToString@x[[2]],{"+","-"}]:>
	{StringJoin@@x,"1"};

(* Optional: Remove charge of ions *)
If[OptionValue["IgnoreIonCharge"],
elements=MapAt[
StringDelete[#,{DigitCharacter,"+","-"}]&,
elements,{All,1}]
];

(* Confirm that tally numbers are expressions *)
elements=MapAt[ToExpression,elements,{All,2}];

(* Optional: Keep tally of the various atoms *)
If[!tallyQ,elements=DeleteDuplicates[First/@elements]];

elements
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
GetLatticeParameters::input="Invalid input.";
GetLatticeParameters::space="\"Space\" must either be \"Direct\", \"Reciprocal\" or \"Both\".";

Options@GetLatticeParameters={
"RoundAnglesThreshold"->0.001,
"Space"->"Direct",
"Units"->False
};

SyntaxInformation@GetLatticeParameters={
"ArgumentsPattern"->{_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
GetLatticeParameters[input_,OptionsPattern[]]:=Block[{
ExtractParametersFromMatrix,RoundAngles,UseUnits,
space=OptionValue["Space"],
a,b,c,\[Alpha],\[Beta],\[Gamma],cellDirect,cellReciprocal,cell,
\[Delta]=OptionValue["RoundAnglesThreshold"],fr,
temp},

(* Auxiliary functions *)
ExtractParametersFromMatrix[matrix_]:=(
{a,b,c}=Sqrt@Diagonal@matrix;
\[Alpha]=ArcCos[matrix[[2,3]]/(b*c)]/Degree;
\[Beta]=ArcCos[matrix[[1,3]]/(a*c)]/Degree;
\[Gamma]=ArcCos[matrix[[1,2]]/(a*b)]/Degree;
Return[{a,b,c,\[Alpha],\[Beta],\[Gamma]}]
);

RoundAngles[cell_List]:=(temp=cell;Do[
fr=FractionalPart@temp[[i]];
If[fr>0.5,fr=1-fr];
If[fr<=\[Delta],temp[[i]]=Round@temp[[i]]],
{i,4,6}];
Return@temp);

UseUnits[cell_]:=(temp=cell;Do[
Which[
i<=3,temp[[i]]=Quantity[temp[[i]],"Angstroms"],
i>=4,temp[[i]]=Quantity[temp[[i]],"Degrees"]],
{i,6}];
Return@temp);

(* Optional: Reciprocal space counterparts *)
If[!MemberQ[{"Direct","Reciprocal","Both"},space],
Message[GetLatticeParameters::space];
Abort[]
];

(* Obtain lattice parameters {a,b,c,\[Alpha],\[Beta],\[Gamma]} *)
Which[
(* A. Crystal entry input *)
StringQ[input],	
InputCheck["CrystalQ",input];
cellDirect=QuantityMagnitude@Values@$CrystalData[
input,"LatticeParameters"];
If[space==="Reciprocal"||space==="Both",
cellReciprocal=ExtractParametersFromMatrix@Inverse@GetCrystalMetric@cellDirect
],

(* B. Metric input *)
MatrixQ[input]&&Dimensions[input]=={3,3},	
Which[
space==="Direct",
cellDirect=ExtractParametersFromMatrix[input],

space==="Reciprocal",
cellReciprocal=ExtractParametersFromMatrix[Inverse@input],

space==="Both",
cellDirect=ExtractParametersFromMatrix[input];
cellReciprocal=ExtractParametersFromMatrix[Inverse@input]
],

(* C. None of the above *)
True,	
Message[GetLatticeParameters::input];
Abort[]
];

(* Check options *)
If[space==="Both",
cell={cellDirect,cellReciprocal};
cell=RoundAngles/@cell;
If[TrueQ@OptionValue["Units"],cell=UseUnits/@cell],

cell=If[space==="Direct",cellDirect,cellReciprocal];
cell=RoundAngles@cell;
If[TrueQ@OptionValue["Units"],cell=UseUnits@cell]
];

cell
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
GetLaueClass::missing="No data found on \[LeftGuillemet]`1`\[RightGuillemet].";

SyntaxInformation@GetLaueClass={
"ArgumentsPattern"->{_}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
GetLaueClass[symbol_String]:=Block[{g,extract},
(* $CrystalData *)
g=$CrystalData[symbol]["SpaceGroup"];
If[StringQ@g,Goto["End"]];

(* Special case: -3m *)
Which[
symbol=="\!\(\*OverscriptBox[\(3\), \(_\)]\)1m"||symbol=="-31m",Return["\!\(\*OverscriptBox[\(3\), \(_\)]\)1m"],
symbol=="\!\(\*OverscriptBox[\(3\), \(_\)]\)m:r"||symbol=="-3m:r",Return["\!\(\*OverscriptBox[\(3\), \(_\)]\)m:r"],
True,Null
];

(* Point group or space group *)
g=$GroupSymbolRedirect[symbol]["LaueClass"];
If[StringQ@g,Return@g];

(* Point group or space group (alternative setting) *)
extract=FullForm[Quiet@Extract[
$GroupSymbolRedirect,
symbol,Inactivate]][[1,1]];
g=$GroupSymbolRedirect[extract]["LaueClass"];

(* No match *)
If[!StringQ@g,
Message[GetLaueClass::missing,symbol];Abort[],
Return@g];

(* End (for use with $CrystalData) *)
Label["End"];
GetLaueClass[g]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
GetScatteringCrossSections::"invsrc"="Invalid source.";
GetScatteringCrossSections::"invelt"="Invalid element.";
GetScatteringCrossSections::invproc="The scattering process type \[LeftGuillemet]`1`\[RightGuillemet] is not recognised.";
GetScatteringCrossSections::"invwlrange"="The wavelength, `1` \[CapitalARing], must be within (0.001 \[LessEqual] \[Lambda] \[LessEqual] 3.000) \[CapitalARing] when using cross sections.";

Options@GetScatteringCrossSections={
"PhysicalProcess"->"",
"Source"->"xraylib",
"Units"->True
};

SyntaxInformation@GetScatteringCrossSections={
"ArgumentsPattern"->{_,_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
GetScatteringCrossSections[input_,
wavelength:_?(NumericQ[#]||QuantityQ[#]&):-1,OptionsPattern[]]:=Block[{
src,unitsQ,elements,\[Lambda]=wavelength,pp=OptionValue["PhysicalProcess"],column,
setpos,\[Sigma],file,read,\[Sigma]i},

(*---* Checks *---*)
(* Chemical element(s) *)
elements=Flatten@GetElements[input];

(* Data source *)
src=FileNameJoin[{$MaXrdPath,"Core","Data","CrossSections",OptionValue["Source"]}];
If[!DirectoryQ@src,Message[GetScatteringCrossSections::"invsrc"];Abort[]];
unitsQ=OptionValue["Units"];

(* Wavelength and its range *)
If[!(0.001<=\[Lambda]<=3.000),
Message[GetScatteringCrossSections::invwlrange,ToString@\[Lambda]];
Abort[]];
\[Lambda]=InputCheck["ProcessWavelength","",wavelength];

(* Column to read from (cross section type) *)
column=Which[
pp==="",5,
MemberQ[{
"Photoelectric","Photoionisation"},pp],2,
MemberQ[{
"Coherent","Rayleigh","Thompson",
"Classical","Elastic"},pp],3,
MemberQ[{
"Incoherent","Compton","Inelastic"
},pp],4,
pp==="Total",5,
True,Message[GetScatteringCrossSections::invproc,pp];Abort[]
];

(*---* Read from file *---*)
(* Stream position in dat files *)
setpos=66*(Round[1000*\[Lambda]]-1);

(* Extract cross sections; \[Sigma] = \[Sigma](element, wavelength) *)
\[Sigma]={};
Do[
file=OpenRead@FileNameJoin[{src,X<>".dat"}];
SetStreamPosition[file,setpos];
read=Read[file,Record];
\[Sigma]i=ToExpression@StringReplace[StringSplit[read][[column]],
m__~~"E"~~s:{"+","-"}~~e:DigitCharacter..:>
m<>"*10^("<>s<>e<>")"];
If[unitsQ,\[Sigma]i=Quantity[\[Sigma]i,"Barns"]];
Close[file];
AppendTo[\[Sigma],X->\[Sigma]i],
{X,elements}];

(*---* Output *---*)
Association@\[Sigma]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
GetSymmetryData::InvalidLabel="\[LeftGuillemet]`1`\[RightGuillemet] is not a recognised label.";
GetSymmetryData::Incompatible="Incompatible group type and label.";

Options@GetSymmetryData={
"UnambiguousSymbol"->True,
"UseMainEntry"->False
};

SyntaxInformation@GetSymmetryData={
"ArgumentsPattern"->{_,_.,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
GetSymmetryData[input_String,label_String:"Lookup",OptionsPattern[]]:=Block[{
group,validLabels,type,data,dataMain,temp},

(* Extract point- or space group (also check $CrystalData) *)
group=InputCheck["GetPointSpaceGroupCrystal",input];
validLabels={"Lookup","Symbol","HermannMauguinFull","HermannMauguinShort",	"HallString","PointGroupNumber","SpaceGroupNumber","LaueClass","CrystalSystem","Centring","MainEntryQ","GroupType","Setting"};

If[!MemberQ[validLabels,label],
Message[GetSymmetryData::InvalidLabel,label];Abort[]];

data=$GroupSymbolRedirect@group;
type=If[MemberQ[$PointGroups,data,Infinity],"PointGroup","SpaceGroup"];

If[
(type==="PointGroup"&&label==="SpaceGroupNumber")||
(type==="SpaceGroup"&&label==="PointGroupNumber"),
Message[GetSymmetryData::incompatible];Abort[]];

If[label==="GroupType",Return@type];

(* Group designation *)
group=If[type==="SpaceGroup",
data["Name","HermannMauguinFull"],
data["Name","Symbol"]
];

If[label=="Centring",
If[type=="PointGroup",Message[GetSymmetryData::Incompatible];Abort[]];
Return@StringTake[group,1]
];

(* Ascertain main entry *)
dataMain=If[type==="SpaceGroup",	
temp=Position[$SpaceGroups,group];
$SpaceGroups[temp[[1,1,1]]],

temp=Position[$PointGroups,data];
$PointGroups[temp[[1,1,1]]]
];

(* Optional: Use main entry corresponding to input *)
If[TrueQ@OptionValue["UseMainEntry"],data=dataMain];

(* Executing commands *)
If[label==="Lookup",Return@data];

If[label==="MainEntryQ",Return@KeyExistsQ[data,type<>"Number"]];

If[MemberQ[{
"Symbol","HermannMauguinFull",
"HermannMauguinShort","HallString"},label],
	
(* Optional: Let output be unambiguous *)
If[TrueQ@OptionValue["UnambiguousSymbol"]&&label==="Symbol",
Return@ToStandardSetting@data["Name","HermannMauguinFull"],
Return@data["Name",label]]
];

If[label==="Setting",Return@data@label];

If[MemberQ[{
"PointGroupNumber","SpaceGroupNumber",
"LaueClass","CrystalSystem"},label],
Return@dataMain@label
]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
GetSymmetryOperations::missing="No data found on \[LeftGuillemet]`1`\[RightGuillemet].";

Options@GetSymmetryOperations={
"UseCentring"->False
};

SyntaxInformation@GetSymmetryOperations={
"ArgumentsPattern"->{_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
GetSymmetryOperations[input_String,OptionsPattern[]]:=
Block[{temp1,temp2,c},
(* Input check *)
temp1=$GroupSymbolRedirect@InputCheck[
"GetPointSpaceGroupCrystal",input];

(* Point group, alternative setting *)
If[KeyExistsQ[temp1,"MatrixOperations"],
Return@temp1["MatrixOperations"]];

temp2=temp1["SymmetryOperations"];

(*---* Point group *---*)
If[KeyExistsQ[temp2,"MatrixOperations"],
Return@temp2["MatrixOperations"]];

(*--*- Space group *---*)
If[OptionValue["UseCentring"],
(* Apply centring vectors *)
c=StringTake[temp1["Name","HermannMauguinShort"],1];
c=InputCheck["GetCentringVectors",c];
temp2=Table[{
temp2[[i,1]],
Mod[temp2[[i,2]]+c[[j]],1]
},{j,Length@c},{i,Length@temp2}];
Flatten[temp2,1],

(* No centring applied *)
temp2]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
ImportCrystalData::subdataInteger="\"\!\(\*
StyleBox[\"ExtractSubdata\", \"Program\"]\)\" must be a positive integer.";
ImportCrystalData::subdataLength="The \!\(\*
StyleBox[\".\", \"Program\"]\)\!\(\*
StyleBox[\"cif\", \"Program\"]\) file has a subdata length of `1`.";
ImportCrystalData::latticeParameters="No lattice parameters were located, or they have an invalid form.";
ImportCrystalData::atomData="No atom data was located.";
ImportCrystalData::SG="Could not determine space group. 'P1' will be used.";
ImportCrystalData::cell="Could not work out the unit cell properly.";
ImportCrystalData::notMaXrd="Data collected using `1` radiation. Errors may occur.";
ImportCrystalData::modulation="Modulated structure detected. Errors may occur.";

Options@ImportCrystalData={
"DataFile"->FileNameJoin[{$MaXrdPath,"UserData","CrystalData.m"}],
"ExtractSubdata"->1,
"IgnoreIonCharge"->False,
"Notes"-><||>,
"RoundAnglesThreshold"->0.001,
"Units"->True,
"OverwriteWarning"->True
};

SyntaxInformation@ImportCrystalData={
"ArgumentsPattern"->{___,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
ImportCrystalData[
{CrystalName_String,
ChemicalFormula:_?StringQ:"",
Z:_?IntegerQ:0,
SpaceGroup_:1,
Wavelength:_?NumericQ:-1
},
GetLatticeParameters_List,
AtomData_List,
OptionsPattern[]
]:=Block[
{choice,name,sg,cell,\[Delta],fr,latticeItem,\[Lambda],itemAtomData,item,dataFile=OptionValue["DataFile"],temp},

(*---* Check if name already exists *---*)
If[CrystalName==="",
name="ImportedCrystal_"<>DateString["ISODate"],name=CrystalName];
If[OptionValue["OverwriteWarning"],
	If[KeyExistsQ[$CrystalData,name],
	choice=ChoiceDialog["\[LeftGuillemet]"<>name<>
"\[RightGuillemet] already exists in $CrystalData.\nDo you want to overwrite this entry?"]];
	If[!choice,Return[]]
];

(*--- Space Group *---*)
sg=InputCheck["GetPointSpaceGroupCrystal",SpaceGroup];
sg=ToStandardSetting@sg;

(*---* Lattice parameters *---*)
If[!AllTrue[GetLatticeParameters,NumericQ[#]&]||
Length@GetLatticeParameters!=6,
Message[ImportCrystalData::latticeParameters]];

	(* Optional: Round angles *)
	cell=GetLatticeParameters;
	\[Delta]=OptionValue["RoundAnglesThreshold"];
	Do[
	fr=FractionalPart@cell[[i]];
	If[fr>0.5,fr=1-fr];
	If[fr<=\[Delta],cell[[i]]=Round@cell[[i]]],
	{i,4,6}];

	(* Optional: Use units *)
	If[OptionValue["Units"],
	Do[
	Which[
	i<=3,cell[[i]]=Quantity[cell[[i]],"Angstroms"],
	i>=4,cell[[i]]=Quantity[cell[[i]],"Degrees"]],
	{i,6}]];

	(* Prepare association entry *)
	latticeItem=Association@
	Thread[{"a","b","c","\[Alpha]","\[Beta]","\[Gamma]"}->cell];

(*---* Wavelength *---*)
	(* Optional: Use units *)
	If[OptionValue["Units"],
	\[Lambda]=Quantity[Wavelength,"Angstroms"],
	\[Lambda]=Wavelength];

(*---* Atom data *---*)
If[AtomData==={},
itemAtomData={<||>},

itemAtomData=Table[
DeleteMissing@Part[AtomData[[i]],{
"Element",
"OccupationFactor",
"SiteSymmetryMultiplicity",
"SiteSymmetryOrder",
"FractionalCoordinates",
"DisplacementParameters",
"Type"}],
{i,Length@AtomData}];

Do[
If[KeyExistsQ[First@itemAtomData,k],
itemAtomData=MapAt[
ToExpression,itemAtomData,{All,Key[k]}]],
	{k,{
	"OccupationFactor",
	"SiteSymmetryMultiplicity",
	"SiteSymmetryOrder",
	"DisplacementParameters"}}
];

	(* Checking strings in coordinates *)
	Do[itemAtomData[[i,"FractionalCoordinates"]]=itemAtomData[[i,"FractionalCoordinates"]]/.x_String:>ToExpression[x],
	{i,Length@itemAtomData}]
];

(*---* Preparing item *---*)
item=<|
	"ChemicalFormula"->ChemicalFormula,
	"FormulaUnits"->Z,
	"SpaceGroup"->sg,
	"LatticeParameters"->latticeItem,
	"Wavelength"->\[Lambda],
	"AtomData"->itemAtomData,
	"Notes"->OptionValue["Notes"]
	|>;

(* Delete certain keys *)
	If[item["ChemicalFormula"]==="",
	KeyDropFrom[item,"ChemicalFormula"]];

	If[item["Notes"]===<||>,
	KeyDropFrom[item,"Notes"]];

	If[!Positive@item["Wavelength"],
	KeyDropFrom[item,"Wavelength"]];

	If[Z===0,
	KeyDropFrom[item,"FormulaUnits"]];

	(* If all occupation factors = 1, delete column *)
	temp=item[["AtomData",All,"OccupationFactor"]];
	If[AllTrue[N@temp,#==1.&],
	item[["AtomData",All]]=KeyDrop[
	item[["AtomData",All]],"OccupationFactor"]];

	(* If all displacement parameters = 0, delete column *)
	temp=item[["AtomData",All,"DisplacementParameters"]];
	If[AllTrue[N@temp,#==0.&],
	item[["AtomData",All]]=KeyDrop[
	item[["AtomData",All]],
	{"DisplacementParameters","Type"}]];


InputCheck["Update$CrystalDataFile",dataFile,name,item];

(* Display *)
KeyValueMap[
If[#1=="AtomData",
"AtomData"->Shallow[#2,1],
#1->#2]&,
$CrystalData[name]]
]


(* ::Input::Initialization:: *)
ImportCrystalData[ciffile_,Name_String:"",OptionsPattern[]]:=Block[{
(* A. Input check and setup *)name,import,sub,endstring,enc,left,mid,right,
modulationQ=False,
(* B. Lattice parameters *)cell,x,X,multipleQ,parts,coordCount,
(* C. Atom data *)atomdata,atomtags,c,
(* D. Anisotropic displacement parameters (ADPs) *)anisodata,anisoOrder,P,atomoverview,tags,labels,disp,item,
(* E. Misc data labels (wavelength, formula units) *)\[Lambda],Z,
(* F. Chemical formula *)formula,chemicalformula,L,l,r,checkParentheses,
(* G. Space group *)sgTags,sgData,sg,
(* H. Adding item to dataset *)options,
(* Misc *)temp},

(*---* A. Input check and setup *---*)
If[Name=="",name=FileBaseName[ciffile],name=Name];

(* A.1. Check file *)
import=Check[Import[ciffile,"String"],Abort[]];
sub=OptionValue["ExtractSubdata"];
	If[!(IntegerQ[#]&&Positive[#]&@sub),
	Message[ImportCrystalData::subdataInteger];
	Abort[]];

(* A.2. Auxiliary variables *)
endstring={"loop_","\n\n",";","#",EndOfString};
enc={"'","\""};(* annotation/enclosing marks *)
{left,mid,right}={"\!\(\*SubscriptBox[\(","\), \(","\)]\)"};

(* A.3. Check radiation type *)
temp=StringCases[import,
Shortest[{
"_diffrn_radiation_type","_diffrn_radiation_probe"}
~~Whitespace~~{"",enc}~~t:LetterCharacter..~~"\n"]:>t];
If[MemberQ[temp,#],Message[ImportCrystalData::notMaXrd,#]]
&/@{"neutron","electron"};

(* A.4. Check for modulation *)
temp=If[StringContainsQ[import,"_space_group_ssg_name"],
modulationQ=True;
Message[ImportCrystalData::modulation]];

(*---* B. Lattice parameters *---*)
(* B.1. Extracting lattice parameters *)
cell=StringCases[import,Shortest[
x:("_cell_"~~{"length","angle"}~~__~~
{DigitCharacter,"."}..)~~
{"(",Whitespace}]:>ToLowerCase@x,
IgnoreCase->True];

	(* Check *)
	If[cell=={},
	Message[ImportCrystalData::latticeParameters];
	Abort[]];

(* B.2. Check for multiplue structures *)
	cell=StringSplit[cell,Whitespace];
	Which[
	Length@cell>6,
		multipleQ=True;parts=Length[cell]/6,
	Length@cell==6,
		multipleQ=False;parts=1,
	True,
		Message[ImportCrystalData::cell];
		Return@cell
	];

	(* Correct ordering *)
	cell=Partition[cell,
	Length@cell/Quotient[Length@cell,6]];

	Do[
	X=cell[[i]];
	x=X[[All,1]];
	P=FindPermutation[x,{
	"_cell_length_a",
	"_cell_length_b",
	"_cell_length_c",
	"_cell_angle_alpha",
	"_cell_angle_beta",
	"_cell_angle_gamma"}];
	cell[[i]]=Permute[X,P],
	{i,Length@cell}];

	cell=cell[[All,All,2]];

(* B.3. Check subdata extraction *)
	(* Verify with fractional coordinates *)
	coordCount=StringCount[import,"_atom_site_fract_x"];
	If[coordCount===0,
	Message[ImportCrystalData::atomData];Abort[]];	
	
	parts=Min[parts,coordCount];

	If[(multipleQ&&sub>parts)||(!multipleQ&&sub!=1),
	Message[ImportCrystalData::subdataLength,parts];
	Abort[]];

	If[!IntegerQ[parts],
	Message[ImportCrystalData::subdataInteger];
	Return@parts];

	(* Extract subdata *)
	cell=ToExpression@cell[[sub]];


(*---* C. Atom data *---*)
Label["AtomData"];

(* C.1. Extracting relevant data block *)
	(* Fractional coordinates *)
	(* Occupation factor *)
	(* Site symmetry multiplicity *)

	(* Extracting _atom_site loop *)
	atomdata=StringCases[import,Shortest[
labels:(Whitespace~~"_atom_site_"~~__~~"\n")~~data:(StartOfLine~~Whitespace...~~LetterCharacter~~__~~"\n")
~~{endstring,"_atom_site_aniso","_"~~Except["a"]}]:>{labels,data}];

		(* No data? *)
		If[atomdata=={},
		Message[ImportCrystalData::atomData];
		Abort[]];

		(* Delete cases containing anisotropy data *)
		atomdata=DeleteCases[atomdata,
		x_/;StringContainsQ[x[[1]],"aniso"]];

	(* Specify sub-data *)
	atomdata=atomdata[[sub]];

(* C.2. Organising data *)
atomtags=Flatten@StringCases[atomdata[[1]],
"_atom_site_"~~{WordCharacter,"_"}..];

atomdata=StringDelete[atomdata[[2]],"("~~DigitCharacter..~~")"];
atomdata=Partition[StringSplit@atomdata,Length@atomtags];
atomdata=DeleteCases[atomdata,x_/;Length[x]!=Length[atomtags]];
atomdata=Association@Thread[atomtags->Transpose@atomdata];

(* C.3 Fixing entries *)
(* If 'site_type_symbol' is missing, copy 'site_label' *)
If[!KeyExistsQ[atomdata,"_atom_site_type_symbol"],
AppendTo[atomdata,
"_atom_site_type_symbol"->atomdata["_atom_site_label"]]];

(* Process and check elements *)
temp=atomdata["_atom_site_type_symbol"];
temp=InputCheck["InterpretElement",temp];

(* Optional: Clear any ion charges *)
If[OptionValue["IgnoreIonCharge"],
temp=StringDelete[temp,{"+","-",DigitCharacter}]];

(* Update 'atomdata' with 'temp' *)
atomdata["_atom_site_type_symbol"]=temp;


(*---* D. Anisotropic displacement parameters *---*)
(* D.1. If missing, use default values for ADP *)
L=Length@First@atomdata;
If[!KeyExistsQ[atomdata,"_atom_site_adp_type"],
AppendTo[atomdata,"_atom_site_adp_type"->ConstantArray["Uiso",L]]];
If[!KeyExistsQ[atomdata,"_atom_site_U_iso_or_equiv"],
AppendTo[atomdata,"_atom_site_U_iso_or_equiv"->ConstantArray[0,L]]];

(* D.2. Anisotropic displacement parameters *)
anisodata=StringCases[import,
Shortest["loop_"~~Whitespace~~
x:("_atom_site_aniso"~~__)~~endstring]:>x];

	(* Check *)
	If[anisodata==={},
	Goto["OrganiseAtomdata"],
	anisodata=anisodata[[sub]]];

(* D.3. Noting the order (permutation) *)
anisoOrder=Flatten@StringCases[anisodata,"U_"~~DigitCharacter..];
P=FindPermutation[{"U_11","U_22","U_33","U_12","U_13","U_23"},anisoOrder];

	(* Nothing there to extract? *)
	temp=Flatten@StringCases[anisodata,
	Shortest["_atom_site_aniso_"~~x__~~Whitespace]:>x];
	c=Flatten@Quiet@Position[temp,
	x_/;StringContainsQ[x,"U"]];

(* D.4. Extracting relevant data and trimming *)
anisodata=StringCases[anisodata,
Shortest["_atom_site_aniso"~~__~~EndOfLine]~~Whitespace~~(* Last line *)
x:(WordCharacter~~__):>x(* Content *)];

(* Check if there is any actual data *)
If[anisodata==={},Goto["OrganiseAtomdata"]];
anisodata=StringSplit[First@anisodata,Whitespace];
anisodata=StringReplace[anisodata,x__~~"("~~__~~")":>x];
anisodata=Partition[anisodata,Length@temp];

(* Correcting parameter order *)
anisodata[[All,c]]=Permute[#,P]&/@anisodata[[All,c]];

(* Associating each atom with values *)
anisodata=Association[
Table[anisodata[[i,1]]->anisodata[[i,c]],
{i,Length@anisodata}]];

(* D.5. Organising the atom data *)
Label["OrganiseAtomdata"];
atomoverview={};
tags={
"_atom_site_occupancy",
"_atom_site_site_symmetry_multiplicity",
"_atom_site_site_symmetry_order"};
labels={
"OccupationFactor",
"SiteSymmetryMultiplicity",
"SiteSymmetryOrder"};

Do[
item=<||>;
AppendTo[item,
"Element"->atomdata[["_atom_site_type_symbol",i]]];
Do[
If[KeyExistsQ[atomdata,tags[[j]]],
AppendTo[item,labels[[j]]->
atomdata[[tags[[j]],i]]]],
{j,Length@tags}];

AppendTo[item,"FractionalCoordinates"->
Evaluate[atomdata[["_atom_site_fract_"<>#,i]]&/@{"x","y","z"}]];
If[StringTake[
atomdata[["_atom_site_adp_type",i]],-3]==="ani",
disp=Part[anisodata,atomdata[["_atom_site_label",i]]],
disp=atomdata[["_atom_site_U_iso_or_equiv",i]]];

AppendTo[item,"DisplacementParameters"->disp];
AppendTo[item,"Type"->atomdata[["_atom_site_adp_type",i]]];
AppendTo[atomoverview,item],{i,Length@First@atomdata}];


(*---* E. Misc data labels *---*)
(* E.1. Wavelength *)
\[Lambda]=StringCases[import,Shortest["_diffrn_radiation_wavelength"~~
Whitespace~~x:{DigitCharacter,"."}..~~Whitespace]:>x];

If[\[Lambda]=={},\[Lambda]=-1,
If[Length@\[Lambda]>1,
\[Lambda]=ToExpression@\[Lambda][[sub]],
\[Lambda]=ToExpression@First@\[Lambda]]];

(* E.2. Forumla units (Z) *)
Z=StringCases[import,Shortest["_cell_formula_units_Z"~~
Whitespace~~z:DigitCharacter]:>z];

If[Z==={},Z=0,Z=ToExpression@First@Z];


(*---* F. Chemical formula *---*)
(* F.1. Extracting formula *)
formula=StringCases[import,
Shortest[#~~{Whitespace,"\n"}..~~{"'","\""}~~f__~~
{"'","\""}~~{Whitespace,"\n"}..]:>f]&/@{
(* Prioritised order *)
"_chemical_formula_iupac",
"_chemical_formula_structural",
"_chemical_formula_sum"};

formula=Select[Flatten@formula,!StringContainsQ[#,{",","?"}]&];
formula=StringDelete[formula,"\r"];

(* F.2. Check for simplest formula *)
temp=Select[formula,!StringContainsQ[#,"("]&];
If[temp=!={},formula={First@temp}];

(* F.3. Misc treatmeant and possible subdata selection *)
If[formula==={}||formula==={""},
chemicalformula="";Goto["SpaceGroup"]];

formula={StringDelete[StringTrim@First@formula,{"'","\""}]};

If[formula==={""},
chemicalformula="";Goto["SpaceGroup"],
formula=StringSplit@formula];

If[Length@formula>1,
formula=formula[[sub]],
formula=Flatten@formula];

(* F.4. Loop for formatting the chemical formula *)
Label["FormatFormula"];
If[AnyTrue[formula,StringContainsQ[#,{"(",")"}]&],
{l,r}=Flatten@Position[StringPosition[
formula,{"(",")"}],{{_,_}}];
	checkParentheses=True];

chemicalformula={};
Do[
temp=Flatten@StringCases[formula[[i]],
x:LetterCharacter..~~
y:({DigitCharacter,"."})..:>{x,y}];
Which[
temp=={},AppendTo[chemicalformula,formula[[i]]],
temp[[2]]=="1",AppendTo[chemicalformula,temp[[1]]],
True,AppendTo[chemicalformula,left<>temp[[1]]<>mid<>temp[[2]]<>right]],
{i,Length@formula}];

chemicalformula=StringDelete[chemicalformula,{"(",")"}];

(* Adding back parentheses *)
If[checkParentheses,
chemicalformula[[l]]="("<>chemicalformula[[l]];
chemicalformula[[r]]=chemicalformula[[r]]<>")"
];

chemicalformula=StringJoin@chemicalformula;


(*---* G. Space group *---*)
Label["SpaceGroup"];

(* G.1. Prioritised list of data labels *)
sgTags={
"_space_group_name_Hall",
"_space_group_name_H-M",
"_space_group_IT_number",
"_symmetry_space_group_name_Hall",
"_symmetry_space_group_name_H-M",
"_symmetry_Int_Tables_number"};

(* G.2. Extract space group sections from imported data *)
sgData=StringCases[import,Shortest[sgTags~~__~~"loop_"]];
If[Length[sgData]>0,sgData=sgData[[sub]]];
sgData=StringTrim@StringDelete[sgData,"loop_"];

(* G.3. Make association of tags and corresponding info *)
sgData=StringCases[sgData,t:sgTags~~Whitespace~~sg:Shortest[Except[WhitespaceCharacter]~~__]~~{"\n",EndOfString}:>{t->sg}];
sgData=Association@Flatten@sgData;
sgData=DeleteCases[sgData,"?"];

(* G.4. Go through priority order and validate *)
Do[
sg=sgData[sgTags[[i]]];
sg=Quiet@InputCheck["InterpretSpaceGroup",sg,False];
If[sg=!=Null,Break[]],
{i,Length@sgTags}];

(* G.5. For modulated structures *)
If[sg===Null&&modulationQ,
sg=StringCases[import,
"_space_group_ssg_name"~~Whitespace~~
{"'","\""}~~sg__~~{"'","\""}~~"\n":>sg];
sg=StringCases[sg,Shortest[StartOfString~~___
~~sg:(LetterCharacter~~__)~~"("]:>sg];
If[sg=!={},sg=Quiet@InputCheck["InterpretSpaceGroup",
First@Flatten@sg,False]]
];

(* G.6. If missing space group, display message and use 'P1' *)
If[sg===Null,Message[ImportCrystalData::SG];sg="P1"];


(*---* H. Adding item to dataset *---*)
options=Thread[#->OptionValue[#],String]&/@(
First/@Options@ImportCrystalData);

ImportCrystalData[
{name,chemicalformula,Z,sg,\[Lambda]},
cell,atomoverview,
options]
]


(* ::Input::Initialization:: *)
ImportCrystalData["RunDialogue"]:=DialogInput[DynamicModule[{
name,gridA,gridB,currentGrid,updGrid,
sgKeys=Keys@$SpaceGroups,sgNumber=Null,sgSymbol=Null,
crystalSystem,systemParameterFields,parameterFields,
a=Null,b=Null,c=Null,
\[Alpha]=Null,\[Beta]=Null,\[Gamma]=Null,
chemicalFormula=Null,wavelength=Null,massDensity=Null,formulaUnits=Null,
atomdata={},atomdataSummary,updAtomdataSummary,
createDeleteButtons,deleteButtons,createAtomdataPanel,atomdataPanel="(no entries)",updAtomdataPanel,
(* Adding entries *)
element=Null,elementListWithNumber,
coordinates,
coordX=Null,coordY=Null,coordZ=Null,
occupationFactor=Null,
adpType="Isotropic",createAdpField,adpField,
adpU11=Null,adpU22=Null,adpU33=Null,
adpU12=Null,adpU13=Null,adpU23=Null,
adpUiso=Null,ADPs,
updRet,toBeReturned,validQ=False
},

(*---* Setup for grid B *---*)
(* Logic *)
elementListWithNumber=#<>" ("<>
ToString@$PeriodicTable[#,"AtomicNumber"]<>
")"&/@Keys@$PeriodicTable;
elementListWithNumber=Thread[
Keys@$PeriodicTable->elementListWithNumber];

createAdpField[]:=Which[
adpType==="Isotropic",
adpField=Column[{
InputField[Dynamic[adpUiso,(adpUiso=#;updRet[])&],
Number,FieldSize->{5.,1.},FieldHint->"\!\(\*SubscriptBox[\(U\), \(iso\)]\)"],
Spacer[{165.,17.5}]
}],

adpType==="Anisotropic",
adpField=Column[{
Row[{
InputField[Dynamic[adpU11,(adpU11=#;updRet[])&],
Number,FieldSize->{5.,1.},FieldHint->"\!\(\*SubscriptBox[\(U\), \(11\)]\)"],
Spacer[5],
InputField[Dynamic[adpU22,(adpU22=#;updRet[])&],
Number,FieldSize->{5.,1.},FieldHint->"\!\(\*SubscriptBox[\(U\), \(22\)]\)"],
Spacer[5],
InputField[Dynamic[adpU33,(adpU33=#;updRet[])&],
Number,FieldSize->{5.,1.},FieldHint->"\!\(\*SubscriptBox[\(U\), \(33\)]\)"]
}],
Row[{
InputField[Dynamic[adpU12,(adpU12=#;updRet[])&],
Number,FieldSize->{5.,1.},FieldHint->"\!\(\*SubscriptBox[\(U\), \(12\)]\)"],
Spacer[5],
InputField[Dynamic[adpU13,(adpU13=#;updRet[])&],
Number,FieldSize->{5.,1.},FieldHint->"\!\(\*SubscriptBox[\(U\), \(13\)]\)"],
Spacer[5],
InputField[Dynamic[adpU23,(adpU23=#;updRet[])&],
Number,FieldSize->{5.,1.},FieldHint->"\!\(\*SubscriptBox[\(U\), \(23\)]\)"]
}]
}]
];
createAdpField[];

(* Dynamic variable storing input data *)
ADPs=Dynamic[{adpU11,adpU22,adpU33,adpU12,adpU13,adpU23}];
coordinates=Dynamic[{coordX,coordY,coordZ}];

updRet[]:=(
toBeReturned=Association[];

If[element=!=Null,
AssociateTo[toBeReturned,"Element"->element],
KeyDropFrom[toBeReturned,"Element"]];

If[occupationFactor=!=Null,
AssociateTo[toBeReturned,"OccupationFactor"->occupationFactor],
KeyDropFrom[toBeReturned,"OccupationFactor"]];

If[AllTrue[coordinates[[1]],NumericQ],
AssociateTo[toBeReturned,"FractionalCoordinates"->coordinates],
KeyDropFrom[toBeReturned,"FractionalCoordinates"]];

Which[
(adpType==="Anisotropic")&&AllTrue[ADPs[[1]],NumericQ],
	AssociateTo[toBeReturned,"DisplacementParameters"->ADPs];
	AssociateTo[toBeReturned,"Type"->"Uani"],
(adpType==="Isotropic")&&NumericQ@adpUiso,
	AssociateTo[toBeReturned,"DisplacementParameters"->adpUiso];
	AssociateTo[toBeReturned,"Type"->"Uiso"],
True,
	KeyDropFrom[toBeReturned,"DisplacementParameters"]
];

If[(element=!=Null)&&(AllTrue[coordinates[[1]],NumericQ]),
validQ=True,validQ=False]
);
updRet[];
(*---* // End of grid B setup // *---*)

(* Logic *)
systemParameterFields[system_String]:=Which[
system==="Triclinic",
Column[{
Row[{
InputField[Dynamic@a,Number,
FieldHint->"a",FieldSize->{5.,1.},
Enabled->True],
Spacer[5],
InputField[Dynamic@b,Number,
FieldHint->"b",FieldSize->{5.,1.},
Enabled->True],
Spacer[5],
InputField[Dynamic@c,Number,
FieldHint->"c",FieldSize->{5.,1.},
Enabled->True]
}],
Row[{
InputField[Dynamic@\[Alpha],Number,
FieldHint->"\[Alpha]",FieldSize->{5.,1.},
Enabled->True],
Spacer[5],
InputField[Dynamic@\[Beta],Number,
FieldHint->"\[Beta]",FieldSize->{5.,1.},
Enabled->True],
Spacer[5],
InputField[Dynamic@\[Gamma],Number,
FieldHint->"\[Gamma]",FieldSize->{5.,1.},
Enabled->True]
}]
}],

system==="Monoclinic",
Column[{
Row[{
InputField[Dynamic@a,Number,
FieldHint->"a",FieldSize->{5.,1.},
Enabled->True],
Spacer[5],
InputField[Dynamic@b,Number,
FieldHint->"b",FieldSize->{5.,1.},
Enabled->True],
Spacer[5],
InputField[Dynamic@c,Number,
FieldHint->"c",FieldSize->{5.,1.},
Enabled->True]
}],
Row[{
InputField[\[Alpha]=90,Number,
FieldHint->"\[Alpha]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[\[Beta]=90,Number,
FieldHint->"\[Beta]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[Dynamic@\[Gamma],Number,
FieldHint->"\[Gamma]",FieldSize->{5.,1.},
Enabled->True]
}]
}],

system==="Orthorhombic",
Column[{
Row[{
InputField[Dynamic@a,Number,
FieldHint->"a",FieldSize->{5.,1.},
Enabled->True],
Spacer[5],
InputField[Dynamic@b,Number,
FieldHint->"b",FieldSize->{5.,1.},
Enabled->True],
Spacer[5],
InputField[Dynamic@c,Number,
FieldHint->"c",FieldSize->{5.,1.},
Enabled->True]
}],
Row[{
InputField[\[Alpha]=90,Number,
FieldHint->"\[Alpha]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[\[Beta]=90,Number,
FieldHint->"\[Beta]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[\[Gamma]=90,Number,
FieldHint->"\[Gamma]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]]
}]
}],

system==="Tetragonal",
Dynamic@Column[{
Dynamic@Row[{
InputField[Dynamic[a,(a=#;b=#)&],Number,
FieldHint->"a",FieldSize->{5.,1.},
Enabled->True],
Spacer[5],
InputField[b=Dynamic@a,Number,
FieldHint->"b",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[Dynamic@c,Number,
FieldHint->"c",FieldSize->{5.,1.},
Enabled->True]
}],
Dynamic@Row[{
InputField[\[Alpha]=90,Number,
FieldHint->"\[Alpha]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[\[Beta]=90,Number,
FieldHint->"\[Beta]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[\[Gamma]=90,Number,
FieldHint->"\[Gamma]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]]
}]
}],

system==="Hexagonal"||system==="Trigonal",
Dynamic@Column[{
Dynamic@Row[{
InputField[Dynamic[a,(a=#;b=#)&],Number,
FieldHint->"a",FieldSize->{5.,1.},
Enabled->True],
Spacer[5],
InputField[b=Dynamic@a,Number,
FieldHint->"b",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[Dynamic@c,Number,
FieldHint->"c",FieldSize->{5.,1.},
Enabled->True]
}],
Dynamic@Row[{
InputField[\[Alpha]=90,Number,
FieldHint->"\[Alpha]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[\[Beta]=90,Number,
FieldHint->"\[Beta]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[\[Gamma]=120,Number,
FieldHint->"\[Gamma]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]]
}]
}],

system==="Cubic",
Dynamic@Column[{
Dynamic@Row[{
InputField[Dynamic[a,(a=#;b=#;c=#)&],Number,
FieldHint->"a",FieldSize->{5.,1.},
Enabled->True],
Spacer[5],
InputField[b=Dynamic@a,Number,
FieldHint->"b",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[c=Dynamic@a,Number,
FieldHint->"c",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]]
}],
Dynamic@Row[{
InputField[\[Alpha]=90,Number,
FieldHint->"\[Alpha]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[\[Beta]=90,Number,
FieldHint->"\[Beta]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]],
Spacer[5],
InputField[\[Gamma]=90,Number,
FieldHint->"\[Gamma]",FieldSize->{5.,1.},
Enabled->False,
BaseStyle->GrayLevel[0.65]]
}]
}]
];

(* Atom data *)
updAtomdataSummary[]:=atomdataSummary=atomdata[[
All,{"Element","FractionalCoordinates"}]];

createDeleteButtons[]:=deleteButtons=
Button["\[Times]",atomdataSummary=Delete[atomdataSummary,#];
createDeleteButtons[];createAtomdataPanel[],
Appearance->"Frameless"]&/@Range@Length@atomdataSummary;

createAtomdataPanel[]:=If[deleteButtons==={},
atomdataPanel="(no entries)",
atomdataPanel=Row[{Pane[
Grid[Transpose@Join[Transpose@Values@atomdataSummary,
{deleteButtons}],Alignment->Left],Scrollbars->False]}]];

updAtomdataPanel[]:=If[atomdata==={},
atomdataPanel="(no entries)",

updAtomdataSummary[];
createDeleteButtons[];
createAtomdataPanel[]
];

(* Initial settings *)
parameterFields=systemParameterFields["Triclinic"];
updAtomdataPanel[];

(*---* Grids *---*)
(* Grid switcher *)
updGrid[label_String]:=currentGrid=Which[
label==="A",gridA,
label==="B",gridB];
updGrid["A"];

(* Grid A *)
gridA=Dynamic@{
{"Crystal name",InputField[Dynamic@name,String,FieldHint->"Crystal name or label"]},
{"Space group",Dynamic@Row[{
PopupMenu[Dynamic[sgNumber,(sgNumber=#;
sgSymbol=sgKeys[[sgNumber]];
crystalSystem=$SpaceGroups[sgSymbol,"CrystalSystem"];
parameterFields=systemParameterFields[crystalSystem]
)&],
	Range@230,"Number"],
PopupMenu[Dynamic[sgSymbol,(sgSymbol=#;
sgNumber=$SpaceGroups[sgSymbol,"SpaceGroupNumber"];
crystalSystem=$SpaceGroups[sgSymbol,"CrystalSystem"];
parameterFields=systemParameterFields[crystalSystem]
)&],
	Keys@$SpaceGroups,"Symbol"]}]
},
{Tooltip["Lattice parameters",Column[{"a  b  c","\[Alpha]  \[Beta]  \[Gamma]","\[ARing]ngstr\[ODoubleDot]m and degree"}],
TooltipDelay->0.6],
Dynamic@parameterFields
},
{},
{Tooltip["Chemical formula","e.g. 'C13 H22 Fe N6 S3'",TooltipDelay->0.6],
InputField[Dynamic@chemicalFormula,String,
FieldHint->"e.g. 'C13 H22 Fe N6 S3'"]},
{Tooltip["Wavelength","\[ARing]ngstr\[ODoubleDot]m",TooltipDelay->0.6],
Row[{
InputField[Dynamic@wavelength,Number,
FieldSize->{5.,1.},FieldHint->"\[Lambda]"],
Spacer[5],
PopupMenu[Dynamic@wavelength,{
1.54059->"Cu \!\(\*SubscriptBox[\(K\[Alpha]\), \(1\)]\)",
1.54443->"Cu \!\(\*SubscriptBox[\(K\[Alpha]\), \(2\)]\)",
1.39223->"Cu \!\(\*SubscriptBox[\(K\[Beta]\), \(1\)]\)",
0.70932->"Mo \!\(\*SubscriptBox[\(K\[Alpha]\), \(1\)]\)",
0.71361->"Mo \!\(\*SubscriptBox[\(K\[Alpha]\), \(2\)]\)",
0.63230->"Mo \!\(\*SubscriptBox[\(K\[Beta]\), \(1\)]\)"
},"(Predefined)"]
}]},
{Tooltip["Mass density","g/\!\(\*SuperscriptBox[\(cm\), \(3\)]\)",TooltipDelay->0.6],
Row[{
InputField[Dynamic@massDensity,Number,
FieldSize->{5.,1.},FieldHint->"\[Rho]"],
Spacer[5],
"formula units",
Spacer[5],
InputField[Dynamic@formulaUnits,Number,
FieldSize->{5.,1.},FieldHint->"Z"]
}]},
{},
{"Atom data",Column[{
Button["Add new element",updGrid["B"]],
Dynamic@atomdataPanel
}]
},
{},
{Null,Row[{CancelButton[],DefaultButton[
DialogReturn[
(* Final checks *)
If[!StringQ@name,name="(no name)"];
If[chemicalFormula===Null,chemicalFormula=""];
If[formulaUnits===Null,formulaUnits=0];
If[sgSymbol===Null,sgSymbol=1];
If[wavelength===Null,wavelength=-1];

(* Return *)
MaXrd`Private`$temp={
{name,chemicalFormula,formulaUnits,sgSymbol,wavelength},
{a,b,c,\[Alpha],\[Beta],\[Gamma]},atomdata
};
MaXrd`Private`$temp=Replace[MaXrd`Private`$temp,
x_Dynamic:>x[[1]],-1]
]]}]}
};

(* Grid B *)
gridB=Dynamic@{
{"Element",PopupMenu[Dynamic[element,(element=#;updRet[])&],
elementListWithNumber,"Element"]},
{"Fractional coordinates",Row[{
InputField[Dynamic[coordX,(coordX=#;updRet[])&],
Number,FieldSize->{5.,1.},FieldHint->"x"],
Spacer[5],
InputField[Dynamic[coordY,(coordY=#;updRet[])&],
Number,FieldSize->{5.,1.},FieldHint->"y"],
Spacer[5],
InputField[Dynamic[coordZ,(coordZ=#;updRet[])&],
Number,FieldSize->{5.,1.},FieldHint->"z"]
}]},
{"Occupation factor",InputField[Dynamic[occupationFactor,(occupationFactor=#;updRet[])&],Number,FieldSize->{5.,1.},FieldHint->"1"]},
{Tooltip["ADP type","Anisotropic displacement parameters",TooltipDelay->0.6],
RadioButtonBar[Dynamic[adpType,(adpType=#;createAdpField[];updRet[])&],
{"Isotropic","Anisotropic"},Method->"Active"]},
{Null,Dynamic@adpField},
{},{},{},{},{},{},{},{},
{Null,Row[{
CancelButton[updGrid["A"]],
DefaultButton["Add",
AppendTo[atomdata,
Replace[toBeReturned,x_Dynamic:>x[[1]],-1]];
updAtomdataPanel[];
updGrid["A"],
Enabled->Dynamic@validQ]
}]}
};

(* Grid on display *)
Dynamic@Grid[currentGrid[[1]],
Spacings->{1,0.5},
Alignment->{Left,Center}]
],

(* Dialogue settings *)
WindowTitle->"Add crystal to $CrystalData",
Modal->True,
WindowSize->{310,All}
];


(* ::Input::Initialization:: *)
ImportCrystalData[]:=Block[{name},

MaXrd`Private`$temp=Null;
ImportCrystalData["RunDialogue"];
If[MaXrd`Private`$temp===Null,Return[]];
name=MaXrd`Private`$temp[[1,1]];

(* Execute ImportCrystalData on input data *)
If[ListQ@MaXrd`Private`$temp,
If[!AllTrue[MaXrd`Private`$temp[[2]]/.x_Dynamic:>x[[1]],NumericQ],
Message[ImportCrystalData::latticeParameters];Abort[],
ImportCrystalData@@MaXrd`Private`$temp
]];

(* Reset temporary variable *)
MaXrd`Private`$temp=Null;

(* Display *)
KeyValueMap[
If[#1=="AtomData",
"AtomData"->Shallow[#2,1],
#1->#2]&,
$CrystalData[name]]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
InputCheck::InvalidLabel="\[LeftGuillemet]`1`\[RightGuillemet] is not a recognised check label.";
InputCheck::InvalidTuple="Reflections (and coordinates) must be on a {\!\(\*
StyleBox[\"h\", \"TI\"]\), \!\(\*
StyleBox[\"k\", \"TI\"]\), \!\(\*
StyleBox[\"l\", \"TI\"]\)} (or {\!\(\*
StyleBox[\"x\", \"TI\"]\), \!\(\*
StyleBox[\"y\", \"TI\"]\), \!\(\*
StyleBox[\"z\", \"TI\"]\)}) form";
InputCheck::SingleTupleExpected="Only one `1` expected.";
InputCheck::IntegerExpected="One or more indices are not integers.";
InputCheck::InvalidInputType="Head of indices must be either Integer, String or Symbol.";
InputCheck::MultipleTuplesExpected="At least two reflections are required to make comparisons.";

InputCheck::EnergyUnitExpected="Input does not have a unit compatible with energy or wavelength.";
InputCheck::InvalidEnergyInput="Input must be an energy or wavelength compatible Quantity, or a number.";
InputCheck::EnergyMustBePositive="The wavelength/energy must be positive.";

InputCheck::NotIn$CrystalData="No data found on \[LeftGuillemet]`1`\[RightGuillemet].";
InputCheck::NoWavelengthIncluded="No wavelength was found for crystal \[LeftGuillemet]`1`\[RightGuillemet].";
InputCheck::InvalidUserInput="Invalid user input.";
InputCheck::InvalidPolarisationSetting="Invalid polarisation setting.";

InputCheck::InvalidPointOrSpaceGroup="Unable to interpret \[LeftGuillemet]`1`\[RightGuillemet] as a point- or space group.";
InputCheck::InvalidPointGroup="Unable to interpret \[LeftGuillemet]`1`\[RightGuillemet] as a point group.";
InputCheck::InvalidSpaceGroup="Unable to interpret \[LeftGuillemet]`1`\[RightGuillemet] as a space group.";
InputCheck::InvalidSpaceGroupNumber="Valid space group numbers are between 1 and 230.";
InputCheck::InvalidSpaceGroupInCrystal="Crystal entry \[LeftGuillemet]`1`\[RightGuillemet] has invalid space group \[LeftGuillemet]`2`\[RightGuillemet].";
InputCheck::InvalidSymmetryObject="Unable to interpret \[LeftGuillemet]`1`\[RightGuillemet] as a point group, space group or a crystal.";

InputCheck::InvalidCentring="Invalid space group centring.";
InputCheck::ElementNumber="Element number `1` is out of range.";
InputCheck::ElementFailed="Unable to interpret \[LeftGuillemet]`1`\[RightGuillemet] as a chemical element.";
InputCheck::ElementError="The element \[LeftGuillemet]`1`\[RightGuillemet] cannot be interpreted.";

InputCheck::InvalidCrystalFamily="\[LeftGuillemet]`1`\[RightGuillemet] is not a valid crystal family.";
InputCheck::InvalidDimension="Dimension must be either \"2D\" or \"3D\".";

InputCheck::DomainSizeError="Discrepancy between given domain size and length.";
InputCheck::InvalidRotationPoint="Invalid rotation point.";
InputCheck::InvalidRotationReference="Reference for rotation anchor must either be \"Host\", \"Domain\", \"DomainCentroid\" or \"Unit\".";
InputCheck::InvalidRotationMap="Values of `1`D rotation maps must be `2`.";


SyntaxInformation@InputCheck={
"ArgumentsPattern"->{__}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
InputCheck[input_List,labels___?(SubsetQ[
{"1hkl","1xyz","Integer","Multiple","StringSymbol","WrapSingle"},{#}]&)
]:=Block[{check,hkl,temp},
(* Check labels *)
	check={labels};
	Do[
	temp=check[[i]];
	If[!MemberQ[{
"1hkl","1xyz","Integer","Multiple",
"StringSymbol","WrapSingle"},temp],
	Message[InputCheck::InvalidLabel,temp];Abort[]],
	{i,Length@check}];

(* Dimensions check (always required) *)
	Which[
	(* Single reflection/coordinate *)
	MatchQ[input,{x1_,x2_,x3_}/;!AnyTrue[{x1,x2,x3},ListQ]],
		hkl={input},

	(* Multiple reflections/coordinates *)
	AllTrue[input,
	MatchQ[#,{x1_,x2_,x3_}/;!AnyTrue[{x1,x2,x3},ListQ]]&],
		hkl=input,

	(* None of the above *)
	True,
		Message[InputCheck::InvalidTuple];Abort[]
	];

(* Single reflection/coordinates check *)
	If[MemberQ[check,"1hkl"],
	If[Length@hkl!=1,
	Message[InputCheck::SingleTupleExpected,"reflection"];
	Abort[]]
	];
	If[MemberQ[check,"1xyz"],
	If[Length@hkl!=1,
	Message[InputCheck::SingleTupleExpected,"coordinate"];
	Abort[]]
	];

(* Multiple reflections check *)
	If[MemberQ[check,"Multiple"],
	If[Length@hkl<2,
	Message[InputCheck::MultipleTuplesExpected];
	Abort[]]
	];

(* Integer check *)
	If[MemberQ[check,"Integer"],
	If[!AllTrue[Flatten@hkl,IntegerQ],
	Message[InputCheck::IntegerExpected];
	Abort[]]
	];

(* Check if Integer, String or Symbol *)
	If[MemberQ[check,"StringSymbol"],
	If[!ContainsAll[{Integer,String,Symbol,Times},
	Head/@Flatten@hkl],
	Message[InputCheck::InvalidInputType];
	Abort[]]
	];

(* Wrap single reflections *)
	If[MemberQ[check,"WrapSingle"],
	If[MatchQ[input,
	{_?(!ListQ[#]&),_?(!ListQ[#]&),_?(!ListQ[#]&)}],
	Return[{input}],
	Return[input]
	]
	]
]


(* ::Input::Initialization:: *)
InputCheck["CrystalQ",input_]:=(
(* Check if entry exists in '$CrystalData' *)
If[MissingQ@$CrystalData[input],
Message[InputCheck::NotIn$CrystalData,input];Abort[]]
)


(* ::Input::Initialization:: *)
InputCheck["GetCentringVectors",centring_String]:=Block[{vectors},
Which[
centring==="P",vectors={},
centring==="F",vectors={{1/2,1/2,0},{0,1/2,1/2},{1/2,0,1/2}},
centring==="I",vectors={{1/2,1/2,1/2}},
centring==="R",vectors={{2/3,1/3,1/3},{1/3,2/3,2/3}},
centring==="A",vectors={{0,1/2,1/2}},
centring==="B",vectors={{1/2,0,1/2}},
centring==="C",vectors={{1/2,1/2,0}},
centring==="H",vectors={{2/3,1/3,0},{1/3,2/3,0}},
True,
	Message[InputCheck::InvalidCentring];
	Abort[];
];
PrependTo[vectors,{0,0,0}]
]


(* ::Input::Initialization:: *)
InputCheck["GetCrystalFamilyMetric",family_,dimension_]:=Block[{M,a,b,c,\[Alpha],\[Beta],\[Gamma]},
(* Input checks *)
If[!MemberQ[{
"Cubic","Hexagonal","Tetragonal","Orthorhombic","Monoclinic","Triclinic"},
family],
Message[InputCheck::InvalidCrystalFamily,family];Abort[]];

If[!MemberQ[{"2D","3D"},dimension],
Message[InputCheck::InvalidDimension];Abort[]];

(* Metric *)
M={{1,b Cos[\[Gamma]],c Cos[\[Beta]]},{0,b Sin[\[Gamma]],c (Cos[\[Alpha]]-Cos[\[Beta]] Cos[\[Gamma]]) Csc[\[Gamma]]},{0,0,c *Sqrt[1-Cos[\[Beta]]^2-(Cos[\[Alpha]]-Cos[\[Beta]] Cos[\[Gamma]])^2 Csc[\[Gamma]]^2]}};

M=N[M/.Which[
family==="Cubic",{
a->1.,b->1.,c->1.,
\[Alpha]->90\[Degree],\[Beta]->90\[Degree],\[Gamma]->90\[Degree]
},
family==="Hexagonal",{
a->1.,b->1.,c->1.,
\[Alpha]->90\[Degree],\[Beta]->90\[Degree],\[Gamma]->120\[Degree]
},
family==="Tetragonal",{
a->1.,b->1.,c->1.61803,
\[Alpha]->90\[Degree],\[Beta]->90\[Degree],\[Gamma]->90\[Degree]
},
family==="Orthorhombic",{
a->1.7,b->1.2,c->0.85,
\[Alpha]->90\[Degree],\[Beta]->90\[Degree],\[Gamma]->90\[Degree]
},
family==="Monoclinic",{
a->1.,b->0.7,c->1.2,
\[Alpha]->90\[Degree],\[Beta]->72.\[Degree],\[Gamma]->90\[Degree]
},
family==="Triclinic",{
a->1.3,b->0.8,c->0.9,
\[Alpha]->66.\[Degree],\[Beta]->77.\[Degree],\[Gamma]->88.\[Degree]
}
]];
If[dimension==="2D",M[[{1,2},{1,2}]],M]
]


(* ::Input::Initialization:: *)
InputCheck["GetCrystalFormulaUnits",input_String]:=Block[{output},
(* Check if crystal entry exists *)
InputCheck["CrystalQ",input];
(* Return crystal wavelength if attached *)
If[KeyExistsQ[$CrystalData[input],"FormulaUnits"],
$CrystalData[input,"FormulaUnits"],
(* If not, query user manually *)
output=ToExpression@InputString[
"Cannot determine the number of formula units "<>
"for \[LeftGuillemet]"<>input<>"\[RightGuillemet]."<>"\n"<>
"Please enter that number or the density below."]
];
If[!NumericQ@output,Message[InputCheck::InvalidUserInput];Abort[],
output]
]


(* ::Input::Initialization:: *)
InputCheck["GetCrystalSpaceGroup",input_String]:=Block[{sg},
(* Check if crystal entry exists *)
InputCheck["CrystalQ",input];
(* Check if space group of crystal is valid *)
sg=$CrystalData[input,"SpaceGroup"];
If[!KeyExistsQ[$GroupSymbolRedirect,sg],
Message[InputCheck::InvalidSpaceGroupInCrystal,input,sg];Abort[],
Return@sg]
]


(* ::Input::Initialization:: *)
InputCheck["GetCrystalWavelength",input_String,abortQ_:True]:=(
(* Check if crystal entry exists *)
InputCheck["CrystalQ",input];
(* Return crystal wavelength if attached *)
If[KeyExistsQ[
$CrystalData[input],"Wavelength"],
	Return@$CrystalData[
	input,"Wavelength"],
(* If not, abort OR return '-1' *)
If[abortQ,
Message[InputCheck::NoWavelengthIncluded,input];Abort[],
Return[-1]
]])


(* ::Input::Initialization:: *)
(* Converts to wavelength [\[ARing]ngstr\[ODoubleDot]ms] *)
InputCheck["GetEnergyWavelength",input_,unitsQ_:True]:=Block[{hcKeV=12.398420,\[Lambda]},
(* Only exception *)
If[input===-1,Return[-1]];

(* Check if positive *)
If[!Positive@input,Message[InputCheck::EnergyMustBePositive];Abort[]];

Which[
(* A. Number input *)
NumericQ@input,
	Which[
	(* 1. Assume \[ARing]ngstr\[ODoubleDot]ms *)
	input <= 5.0,\[Lambda]=N@input,
	(* 2. Assume kilo electronvolt *)
	input <= 250.0,\[Lambda]=hcKeV/input,
	(* 3. Assume electronvolts *)
	True,\[Lambda]=1000*hcKeV/input
	],

(* B. Quantity input *)
QuantityQ[input] ,
	(* Convert wavelength or energy to \[ARing]ngstr\[ODoubleDot]ms *)
	Which[
	UnitDimensions[input]==={{"LengthUnit",1}},
		\[Lambda]=UnitConvert[input,"Angstroms"];
		If[unitsQ,Return@\[Lambda],Return@QuantityMagnitude@\[Lambda]],

	CompatibleUnitQ[input,"Joules"],		
		\[Lambda]=hcKeV/QuantityMagnitude@
		UnitConvert[input,"Kiloelectronvols"],

	True,
		Message[InputCheck::EnergyUnitExpected];
		Abort[]
	],

(* C. None of the above *)
True,
	Message[InputCheck::InvalidEnergyInput];Abort[]
];

(* Set in Quantity if desired *)
If[unitsQ,
Quantity[\[Lambda],"Angstroms"],
\[Lambda]]
]


(* ::Input::Initialization:: *)
InputCheck["GetPointSpaceGroupCrystal",input_]:=Block[{
$CrystalDataCombined,sg},
(* Check if space group number is given *)
If[(1<=input<=230)&&IntegerQ@input,Return@$SpaceGroups[[input,"Name","Symbol"]]];
(* If actual input is a point- or space group, return it *)
If[KeyExistsQ[$GroupSymbolRedirect,input],Return@input];
(* Check if crystal entry exists *)
If[AssociationQ@MaXrd`Private`$TempCrystalData,
$CrystalDataCombined=Join[
$CrystalData,MaXrd`Private`$TempCrystalData],
$CrystalDataCombined=$CrystalData];
If[MissingQ@$CrystalDataCombined[input],
Message[InputCheck::InvalidSymmetryObject,input];Abort[]];
(* Check if space group of crystal exists *)
sg=$CrystalDataCombined[input,"SpaceGroup"];
If[!KeyExistsQ[$GroupSymbolRedirect,sg],
Message[InputCheck::InvalidSpaceGroupInCrystal,input,sg];Abort[],
Return@sg]
]


(* ::Input::Initialization:: *)
InputCheck["GetReciprocalImageOrientation",
latticeInput_,hklPlane_,indexLimit_,directSpaceQ_
]:=Block[{
hkl=hklPlane,
abscissaIndex,ordinateIndex,planeConstant,planeIndex,
bottomLeft={-1,-1},bottomRight={1,-1},topRight={1,1},topLeft={-1,1},
M,\[Xi],imageOrientation
},

If[StringQ@hkl,hkl=MillerNotationToList@hkl];

{abscissaIndex,ordinateIndex,planeConstant}={#1,#2,hkl[[#3]]}&@@Flatten[
Position[hkl,#]&/@{"h","k","l",_Integer}];
planeIndex=First@Complement[Range@3,{abscissaIndex,ordinateIndex}];

If[TrueQ@directSpaceQ,
(* a. Direct space vectors (uvw) *)
M=GetCrystalMetric[latticeInput,"ToCartesian"->True];
\[Xi]=2*indexLimit/Max@M;
M=M[[{abscissaIndex,ordinateIndex},{abscissaIndex,ordinateIndex}]];
imageOrientation=\[Xi]*{bottomRight,topLeft,topRight};
imageOrientation=#.M&/@imageOrientation;
imageOrientation=Insert[#,N@planeConstant,planeIndex]&/@imageOrientation;
imageOrientation=Append[#1,#2]&@@@Transpose[
{imageOrientation,{500,500,1}}];
PrependTo[imageOrientation,imageOrientation[[3,{1,2,3}]]],

(* b. Corners in reciprocal space *)
imageOrientation=indexLimit*{bottomLeft,bottomRight,topLeft};
imageOrientation=N@Insert[#,planeConstant,planeIndex]&/@imageOrientation
];

imageOrientation=MapAt[#,imageOrientation,
{All,{abscissaIndex,ordinateIndex}}]&[DecimalForm[
#,{7,4},NumberPadding->{" ","0"}]&];
imageOrientation=Map[ToString,imageOrientation,{2}];
imageOrientation=StringRiffle[#,",  "]&/@imageOrientation
]


(* ::Input::Initialization:: *)
InputCheck["InterpretElement",input_]:=Block[{elementsIn=input,pertiodicTable,elementsRead,elementsReadNeutral,temp},
(*---* A. Input number *---*)
(* A.1. Check whether number is a string *)
If[StringQ@elementsIn,
If[StringMatchQ[elementsIn,NumberString],
elementsIn=ToExpression@elementsIn]];

(* A.2. Check if valid integer (in the periodic table) *)
If[NumericQ@elementsIn,
If[(1<=elementsIn<=Length@$PeriodicTable)&&IntegerQ@elementsIn,Return[(Keys@$PeriodicTable)[[elementsIn]]],
Message[InputCheck::ElementNumber,ToString@elementsIn];Abort[]]
];

(*---* B. Process single string *---*)
(* B.1. Wrap string *)
If[StringQ@elementsIn,elementsIn={elementsIn}];

(*---* C. Process list of strings *---*)
(* C.1. Check if input is a list of strings *)
If[!ListQ@elementsIn,Goto["Failed"]];
If[!AllTrue[elementsIn,StringQ],Goto["Failed"]];

(* C.2. Set of valid symbols from the periodic table *)
pertiodicTable=Keys@$PeriodicTable;

(* C.3. Find (possible) matches, and establish order *)
elementsRead=StringCases[elementsIn,
StartOfString~~a:LetterCharacter~~{"",b:LetterCharacter}~~
{"",n1:DigitCharacter...~~pm:{"+","-"}~~n2:DigitCharacter...}
:>ToUpperCase[a]<>ToLowerCase[b]<>
If[(n1==="")&&(n2===""),"1",""]<>n1<>n2<>pm];

	(* Check for non-elements *)
	If[MemberQ[elementsRead,{}],Message[InputCheck::ElementError,
	Part[elementsIn,Position[elementsRead,{}][[1,1]]]];Abort[]];

	(* Separate element symbol and charge (if any) *)
	elementsRead=Flatten[StringCases[Flatten@elementsRead,
	StartOfString~~a:LetterCharacter..~~b___
	:>{a,If[b==="","&",b]}],1](* '&' should be deleted later *);
	elementsReadNeutral=elementsRead[[All,1]];

(* C.4. Remove second character if not applicable *)
elementsReadNeutral=elementsReadNeutral/.s_String/;
!MemberQ[pertiodicTable,s]:>StringTake[s,1];

	(* Special case: deuterium *)
	elementsReadNeutral=elementsReadNeutral/."D"->"H";
	
(* C.5. Final validity check *)
temp=Complement[elementsReadNeutral,pertiodicTable];
If[temp=!={},Message[InputCheck::ElementError,
Part[elementsIn,Position[elementsReadNeutral,
First@temp][[1,1]]]];Abort[]];

(* C.6. Concatenate elements and charge *)
elementsRead[[All,1]]=elementsReadNeutral;
elementsRead=elementsRead/."1"->"";
elementsRead=StringDelete[StringJoin/@elementsRead,"&"];

Goto["Done"];

(*---* D. Post process *---*)
(* a. Unable to determine chemical element *)
Label["Failed"];
Message[InputCheck::"ElementFailed",ToString@input];
Abort[];

(* b. Return string (or list of strings) *)
Label["Done"];
If[StringQ@input,
elementsRead[[1]],elementsRead]
]


(* ::Input::Initialization:: *)
InputCheck["InterpretSpaceGroup",input_,abortQ_:True]:=Block[{sg=input,o,temp},
(*---* A. Input number *---*)
(* A.1. Check whether number is a string *)
If[StringQ@sg,
If[StringMatchQ[StringTrim@sg,NumberString],
sg=ToExpression@sg]];

(* A.2. Check if valid integer (a canonical space group number) *)
If[NumericQ@sg,
If[(1<=sg<=230)&&IntegerQ@sg,Return@$SpaceGroups[[sg,"Name","Symbol"]],
Message[InputCheck::InvalidSpaceGroupNumber];
If[abortQ,Abort[],Return@Null]]
];

(*---* B. Process string *---*)
(* B.1. Check if string *)
If[!StringQ@sg,Goto["Failed"]];
sg=StringTrim@sg;

(* B.2 Process any annotations *)
If[StringContainsQ[sg,"origin",IgnoreCase->True],
o=ChoiceDialog["Information on cell origin detected."<>
"\n"<>"Please confirm the cell origin.",{1,2},
WindowTitle->"Cell origin"];
sg=StringDelete[sg,Whitespace~~"("~~__~~")"];
];

(* B.3 Tidy string *)
sg=Fold[StringReplace[#1,#2]&,sg,{
(* B.4.1. Remove any _enclosing_ quotation marks *)
(* Note: Some Hall symbols contain double quotation marks *)
StartOfString~~{"'","\""}~~main__~~{"'","\""}~~EndOfString:>main,

(* B.4.2 Uppercase centring 1 *)
StartOfString~~first:{"-"~~_,_}~~rest__~~EndOfString:>
ToUpperCase@first~~ToLowerCase@rest,

(* B.4.3. Fix boxes *)
{"overscriptbox"->"OverscriptBox",
"subscriptbox"->"SubscriptBox"},

(* B.4.4 Uppercase centring 2 *)
"Box[\("~~c_:>"Box[\("<>ToUpperCase[c],

(* B.4.5 Screw axes *)
a:DigitCharacter~~"("~~b:DigitCharacter~~")":>a<>b,

(* B.4.6. Miscellaneous replacements *)
";"->" "
}];
sg=StringTrim@sg;

(*---* C. Check symbol *---*)
(* C.1. Check centring symbol *)
If[!MemberQ[{"P","I","F","R","A","B","C","H","\!","-"},
StringTake[sg,1]],Goto["Failed"]];

(* C.2 Check if found by '$GroupSymbolRedirect' *)
temp=$GroupSymbolRedirect[sg];
If[!MissingQ@temp,
If[KeyExistsQ[temp,"PointGroupNumber"],Goto["Failed"]];
sg=temp[["Name","Symbol"]];Goto["SpaceGroupFound"]];

	(* Exception: Old symbol *)
	If[sg==="Fm3m",sg="Fm-3m";Goto["SpaceGroupFound"]];

(* C.3 Delete whitespace and check again *)
temp=StringDelete[sg,Whitespace];
temp=$GroupSymbolRedirect[temp];
If[!MissingQ@temp,
sg=temp[["Name","Symbol"]];Goto["SpaceGroupFound"]];

(*---* D. Post process *---*)
(* A. Unable to determine space group *)
Label["Failed"];
Message[InputCheck::InvalidSpaceGroup,input];
If[abortQ,Abort[],Return@Null];

(* B. Return non-ambiguous output *)
Label["SpaceGroupFound"];
If[ValueQ[o],sg=sg<>":"<>ToString[o]];

ToStandardSetting[sg]
]


(* ::Input::Initialization:: *)
InputCheck["PointGroupQ",input_String]:=(
(* Check if valid space group string *)
If[!MemberQ[$PointGroups,input,Infinity],
Message[InputCheck::InvalidPointGroup,input];Abort[]
])


(* ::Input::Initialization:: *)
InputCheck["PointSpaceGroupQ",input_String]:=
(* Check if valid point- or space group string *)
	If[!KeyExistsQ[$GroupSymbolRedirect,input],
	Message[InputCheck::InvalidPointOrSpaceGroup,input];Abort[]
	]


(* ::Input::Initialization:: *)
InputCheck["Polarisation",type_String,scatteringAngle_:0]:=
Which[
type==="sigma"||type==="\[Sigma]",1,
type==="pi"||type==="\[Pi]",Abs@Cos[scatteringAngle],
True,Message[InputCheck::InvalidPolarisationSetting];Abort[]
]


(* ::Input::Initialization:: *)
InputCheck["ProcessWavelength",crystal_String,wavelength_,abortQ_:True]:=Block[{\[Lambda]},
If[wavelength===-1,
\[Lambda]=InputCheck["GetCrystalWavelength",crystal,abortQ],
\[Lambda]=wavelength];
InputCheck["GetEnergyWavelength",\[Lambda],False]]


(* ::Input::Initialization:: *)
InputCheck["RotationTransformation",
{{A_Integer,B_Integer,C_Integer},domains_List},
{
AnchorShift_List,
anchorReference_String,
rotations_,
rotationAxes_:IdentityMatrix@3
},
force3Dinterpretation_:False
]:=Module[{
twoDimensionalQ,coordinates,
\[Zeta],R,
uniqueDomains,anchorShift=AnchorShift,anchors,coordinatesGrouped,
domainAnchors,zeroRotation,zeroAnchor
},

(* Preparations and checks *)
If[Length@domains!=A*B*C,
Message[InputCheck::DomainSizeError];Abort[]];

uniqueDomains=DeleteDuplicates@domains;
twoDimensionalQ=(C===1);
If[force3Dinterpretation,twoDimensionalQ=False];
If[twoDimensionalQ,anchorShift=anchorShift[[{1,2}]]];

If[twoDimensionalQ,
If[!MatchQ[anchorShift,{_?NumericQ,_?NumericQ}],
Message[InputCheck::InvalidRotationPoint];Abort[]],
If[!MatchQ[anchorShift,{_?NumericQ,_?NumericQ,_?NumericQ}],
Message[InputCheck::InvalidRotationPoint];Abort[]]
];

If[AssociationQ@rotations,
If[twoDimensionalQ,
If[AnyTrue[Values@rotations,
!NumericQ[#]&],
Message[InputCheck::InvalidRotationMap,"2","scalars"];
Abort[]],
If[AnyTrue[Values@rotations,
!MatchQ[#,ConstantArray[_?NumericQ,3]]&],
Message[InputCheck::InvalidRotationMap,"3","lists of three numbers"];Abort[]]
]
];

If[!MemberQ[
{"Host","Domain","DomainCentroid","Unit"},anchorReference],
Message[InputCheck::InvalidRotationReference];Abort[]];

(* Rotation anchor and reference for domains *)
If[anchorReference==="Domain"||anchorReference==="DomainCentroid",
coordinates=Flatten[Table[{i,j,k},
{i,0.,A-1},{j,0.,B-1},{k,0.,C-1}],2];
If[twoDimensionalQ,coordinates=coordinates[[All,{1,2}]]];
coordinatesGrouped=Pick[coordinates,domains,#]&/@uniqueDomains;
domainAnchors=<|
"Domain"->Flatten[TakeSmallestBy[#,Norm,1]&/@coordinatesGrouped,1],
"DomainCentroid"->(#+If[twoDimensionalQ,{0.5,0.5},{0.5,0.5,0.5}]
&/@Map[Mean,Transpose/@coordinatesGrouped,{2}])
|>@anchorReference;
#+anchorShift&/@domainAnchors;
anchors=N@Association@Thread[uniqueDomains->domainAnchors]
];

(* Rotation transformation function 'R' *)
{zeroRotation,zeroAnchor}=If[twoDimensionalQ,
{0.,{0.,0.}},{{0.,0.,0.},{0.,0.,0.}}];
\[Zeta][d_]:=N@Lookup[rotations,d,zeroRotation];
Which[
anchorReference==="Host"&&twoDimensionalQ,
R[d_]:=Chop@RotationTransform[\[Zeta]@d,anchorShift],

anchorReference==="Host",
R[d_]:=Chop[Composition@@MapThread[RotationTransform[#1,#2,anchorShift]&,
{\[Zeta]@d,rotationAxes}]];
R[d_,p_]:=Chop[Composition@@MapThread[RotationTransform[#1,#2,anchorShift]&,
{If[!AssociationQ@rotations,d,\[Zeta]@d],rotationAxes}]],

anchorReference==="Unit"&&twoDimensionalQ,
R[d_,p_]:=Chop@RotationTransform[\[Zeta]@d,p+anchorShift],

anchorReference==="Unit",
R[d_,p_]:=Chop[Composition@@MapThread[RotationTransform[#1,#2,p+anchorShift]&,
{If[!AssociationQ@rotations,d,\[Zeta]@d],rotationAxes}]],

twoDimensionalQ,
R[d_]:=Chop@RotationTransform[\[Zeta]@d,Lookup[anchors,d,zeroAnchor]+anchorShift],

True,
R[d_,p_]:=Chop[Composition@@MapThread[RotationTransform[#1,#2,
Lookup[anchors,d,zeroAnchor]+anchorShift]&,
{\[Zeta]@d,rotationAxes}]]
];

R
]


(* ::Input::Initialization:: *)
InputCheck["Update$CrystalDataAutoCompletion"]:=(
FE`Evaluate[FEPrivate`AddSpecialArgCompletion[#]]&[
"$CrystalData"->{Keys@$CrystalData,
{"AtomData","ChemicalFormula","FormulaUnits",
"LatticeParameters","Notes","SpaceGroup","Wavelength"}
}];
$CrystalData=KeySort@$CrystalData;
Keys@$CrystalData
)


(* ::Input::Initialization:: *)
InputCheck["Update$CrystalDataFile",
dataFile_String,newStructureLabel_String,hostCopy_]:=(
If[!FileExistsQ@dataFile,$CrystalData=<||>];
AssociateTo[$CrystalData,newStructureLabel->hostCopy];
InputCheck["Update$CrystalDataAutoCompletion"];
Export[dataFile,$CrystalData];
);


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
Options@InterplanarSpacing={
"Units"->True
};

SyntaxInformation@InterplanarSpacing={
"ArgumentsPattern"->{_,_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
InterplanarSpacing[input_,reflections_List,OptionsPattern[]]:=Block[{hkl,H,d,h},

(* Check reflection *)
hkl=InputCheck[reflections,"Integer","WrapSingle"];

(* Reciprocal metric *)
H=GetCrystalMetric[input,"Space"->"Reciprocal"];

(* Interplanar distance *)
d=Reap[Do[
h=hkl[[i]];
Sow[
1/Sqrt[h.H.h]
],
{i,Length@hkl}]
][[2,1]];

(* Option: Units *)
If[OptionValue["Units"],
d=Quantity[d,"Angstroms"]];

(* If only one reflection, return content *)
If[Length@d==1,
First@d,d]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
Options@MergeSymmetryEquivalentReflections={
"ToStandardSetting"->True
};

SyntaxInformation@MergeSymmetryEquivalentReflections={
"ArgumentsPattern"->{_,_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
MergeSymmetryEquivalentReflections[group_String,hkl_List,OptionsPattern[]]:=
Block[{input,sg,merged},

(* Check input *)
	input=InputCheck[hkl,"Integer","WrapSingle"];
	sg=InputCheck["GetPointSpaceGroupCrystal",group];

(* Consider duplicate if they generate same symmetry equivalents *)
	merged=DeleteDuplicatesBy[input,
	Sort@SymmetryEquivalentReflections[sg,#]&];

(* Optional: Use standard setting on indices *)
	If[OptionValue["ToStandardSetting"],
	ToStandardSetting[sg,#]&/@merged,
	merged]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
SyntaxInformation@MillerNotationToList={
"ArgumentsPattern"->{_}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
MillerNotationToList[input_String]:=Block[{L,R,temp,split},
(* Removing overbars *)
L="\!\(\*OverscriptBox[\(";
R="\), \(_\)]\)";
	temp=StringReplace[input,L~~Shortest@x__~~R:>"|"<>"-"<>x<>"|"];

(* Separating indices *)
temp=StringReplace[temp,
(* Sign *)
s:{"","-","|"}~~
{
(* Letters are not joined with digits *)
x:LetterCharacter,
(* Digits could be joined *)
d:DigitCharacter~~
y:Shortest[{
"|",
DigitCharacter..~~{",","|",")"}
}]
}:>"|"<>s<>x<>d<>y<>"|"
];
temp=StringReplace[temp,"|"->","];
temp=StringSplit[temp,","];

	(* Special case: Positive single digits/letters only *)
	If[Length@temp<3,
	temp=StringCases[temp,x:WordCharacter:>x];
	temp=Flatten@DeleteCases[temp,{}]];

(* Trimming *)
temp=StringDelete[temp,{"(",")"}];
temp=DeleteCases[temp,x_/;MemberQ[{"",Null,"{}"},x]];

(* If not three entires, split digits *)
split[x_]:=Flatten@StringCases[x,{
p:{"","-"}~~n:DigitCharacter:>p~~n,
n1:DigitCharacter~~n2:DigitCharacter:>{n1,n2}
}];
temp=temp/.x_List/;Length[x]<3:>split[x];

(* Setting numbers as experssions *)
temp/.x_String/;StringContainsQ[x,DigitCharacter]:>ToExpression[x]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
MillerNotationToString::inv="Invalid index input \[LeftGuillemet]`1`\[RightGuillemet].";

SyntaxInformation@MillerNotationToString={
"ArgumentsPattern"->{{_,_,_}}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
MillerNotationToString[inputRaw_List]:=Block[{L,R,quit,i,H,index,input=inputRaw,presentation,output},
(* Input check *)
Check[InputCheck@inputRaw,Goto["End"]];

(* Shortcuts *)
L="\!\(\*OverscriptBox[\(";
R="\), \(_\)]\)";
quit[index_]:=(
Message[MillerNotationToString::inv,index];
Goto["End"]);

(* Pre-processing input *)
input=input/.x_String/;StringContainsQ[x,"-"]:>
-StringDelete[x,"-"];

(* Converting to string with overbar if negative *)
H={};
Do[
index=input[[i]];
Which[
IntegerQ@index,
	If[index<0,
	AppendTo[H,L<>ToString[-index]<>R],
	AppendTo[H,ToString@index]],
StringQ@index,
	If[StringLength@index!=1,quit[index],
	AppendTo[H,index]],
Head[index]===Real,
	AppendTo[H,ToString@index],
Head[index]===Times,
	If[(index[[1]]===-1)&&(StringQ@index[[2]]),
	If[StringLength@index[[2]]!=1,quit[index],
	AppendTo[H,L<>index[[2]]<>R]],
	quit[index]],
True,
	quit[index]
],{i,3}];

(* Presentation *)
presentation=StringJoin[
"("<>H[[1]]<>"|"<>H[[2]]<>"|"<>H[[3]]<>")"];

(* Only remove commas for single digit integers *)
If[AllTrue[Select[input,NumericQ],(Abs[#]<=9)&&IntegerQ[#]&],
output=StringDelete[presentation,"|"],
output=StringReplace[presentation,"|"->","]];

Return@output;

Label["End"];
input
]


(* ::Input::Initialization:: *)
MillerNotationToString[input_String]:=
MillerNotationToString@MillerNotationToList@input


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
ReciprocalSpaceSimulation::invalid="Invalid input form.";
ReciprocalSpaceSimulation::dep="The layer vectors are not linearly independent.";

Options@ReciprocalSpaceSimulation={
"ReturnData"->False
};

SyntaxInformation@ReciprocalSpaceSimulation={
"ArgumentsPattern"->{_,_.,_,{_,_,_},_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
ReciprocalSpaceSimulation[
crystal_,
lambda:_?(NumericQ[#]||QuantityQ[#]&):-1,
{L1_List,L2_List},origin_List,
res_?(NumericQ[#]&&Positive[#]&),
OptionsPattern[]]:=Block[{
\[Lambda],check,
G,Ginv,B,CrystalDot,CrystalCross,
Hx,Hy,Hz,HCx,HCy,HCz,
U,UB,o,ref,refz,flip,condition,pindex,
\[Epsilon],disk,
hkl,xy,pair,points},

(** Input check **)
	\[Lambda]=InputCheck["ProcessWavelength",crystal,lambda];

	check=Flatten[{L1,L2,origin}];
	If[Length@check!=9||!AllTrue[check,NumericQ],
	Message[ReciprocalSpaceSimulation::invalid];Abort[]];

	(* Check if vectors are linearly independent *)
	If[Det[{{L1.L1,L1.L2},{L2.L1,L2.L2}}]==0,
	Message[ReciprocalSpaceSimulation::dep];Abort[]];

(** Metric information **)
	G=GetCrystalMetric[crystal];
	Ginv=Inverse@G;
	B=CholeskyDecomposition@Inverse@G;

(** Dot- and cross products in reciprocal space **)
	CrystalDot[u_,v_]:=Return[
	Sum[
	Sum[
	Ginv[[i,j]]*u[[i]]*v[[j]],
	{j,3}],
	{i,3}]
	];

	CrystalCross[u_,v_]:=Return[
	Sqrt@Det@Ginv*Table[
	Sum[
	Sum[
	Sum[
	Signature[{i,j,k}]*
	G[[i,t]]*u[[j]]*v[[k]],
	{k,3}],
	{j,3}],
	{i,3}],
	{t,3}]
	];

(** Plane of projection **)
	(* Projection plane in reciprocal space *)
	Hx=L1;
	Hy=L2-Hx*CrystalDot[Hx,L2]/CrystalDot[Hx,Hx];
	Hz=CrystalCross[Hx,Hy];

	(* Components in Cartesian frame *)
	{HCx,HCy,HCz}=Normalize[B.#]&/@{Hx,Hy,Hz};

	(* U and UB matrices for generation of coordinates *)
	U=IdentityMatrix[3].Inverse@Transpose[{HCx,HCy,HCz}];
	UB=Chop[U.B];

	(* Reference position in projection *)
	o=origin;
	ref=UB.o;
	refz=ref[[3]];

	If[Chop[origin]=={0,0,0},
	o=L1+L2;flip=True,
	flip=False];

	condition={
{h_,k_,l_}/;h==#1,
{h_,k_,l_}/;k==#2,
{h_,k_,l_}/;l==#3
}&@@o;
	condition=DeleteCases[condition,c_/;
	If[flip,
	!Equal[c[[2,2]],0],
	   Equal[c[[2,2]],0]
	]
	];

		(* Constant plane index *)
		pindex=condition[[1,2]];
		pindex=ToString@pindex[[1]]->pindex[[2]];

	If[Length@condition>1,
	condition=condition[[All,2]];
	condition={h_,k_,l_}/;#&[And@@condition],
	condition=First@condition];

(** Building **)
	Label["Building"];

	(* Resolution, limits and background *)
	\[Epsilon]=2.0*^-6;
	disk=Show[Graphics[{
	GrayLevel[0.95],
	EdgeForm[Thickness[Small]],
	Disk[{0,0},1/res]
	}],AspectRatio->1];

	(* Generating reflection *)
	hkl=ReflectionList[
	crystal,\[Lambda],condition,
	"HoldIndex"->pindex,
	"SplitEquivalent"->True,
	"ShowProgress"->False];

	(* Generating coordinates *)
	xy=(UB.#-{0,0,refz}&/@hkl)[[All,{1,2}]];
	pair=Transpose[{xy,hkl}];
	pair=Select[pair,
	Sqrt[CrystalDot[#[[2]],#[[2]]]]<1/(1.01*res)&];

	(* Optional: Return data *)
	If[OptionValue["ReturnData"],Return@pair];	

	pair=Tooltip[#1,MillerNotationToString[#2]]
&@@#&/@pair;
	points=ListPlot[pair,
	PlotStyle->PointSize[Large],PlotRange->All];

Show[disk,points,PlotRange->{{-#,#},{-#,#}}&[1.01/res],ImageSize->Large]	
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
ReflectionList::form="Some reflections are not on a {\!\(\*
StyleBox[\"h\", \"TI\"]\), \!\(\*
StyleBox[\"k\", \"TI\"]\), \!\(\*
StyleBox[\"l\", \"TI\"]\)} form.";
ReflectionList::integer="Some reflection indices are not integers.";
ReflectionList::empty="No reflections match the conditions.";
ReflectionList::keep="Invalid setting for the \[LeftGuillemet]Keep\[RightGuillemet] option.";
ReflectionList::index="Invalid index setting.";
ReflectionList::limit="Limit must be a natural number.";

Options@ReflectionList={
"AngleThreshold"->90.*Degree,
"CustomReflections"->False,
"Keep"->All,
"Limit"->30,
"ShowProgress"->True,
"SplitEquivalent"->False,
"HoldIndex"->False,
(* 'ToStandardSetting' options *)
"ToStandardSetting"->True
};

SyntaxInformation@ReflectionList={
"ArgumentsPattern"->{_,_.,_.,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
ReflectionList[
n_?(IntegerQ[#]&&Positive[#]&),
condition___Condition,
OptionsPattern[]]:=Block[
{opt,i,v,h,k,l,hkl},

(* Optional: Hold an index at same value *)
opt=OptionValue["HoldIndex"];
If[TrueQ[Head@opt==Rule],
{i,v}=Part[opt,#]&/@{1,2};
i=StringTake[i,-1],
i="None"];

	(* Check *)
	If[i!="None",
	If[!MemberQ[{"h","k","l"},i]||!IntegerQ[v],
	Message[ReflectionList::index];
	Abort[]]];

hkl=Which[
i=="h",
Flatten[Table[{v,k,l},{k,-n,n},{l,-n,n}],1],

i=="k",
Flatten[Table[{h,v,l},{h,-n,n},{l,-n,n}],1],

i=="l",
Flatten[Table[{h,k,v},{h,-n,n},{k,-n,n}],1],

True,
Tuples[Range[-n,n],3]
];

(* Remove the '000' reflection *)
	hkl=DeleteCases[hkl,{0,0,0}];

(* Sorting *)
	hkl=SortBy[hkl,{Total@Abs[#]&,Negative}];

(* Checking if extra conditions are present *)
	If[{condition}!={},
	(* Filtering reflections *)
	hkl=Cases[hkl,condition]];

Return@hkl
]


(* ::Input::Initialization:: *)
ReflectionList[
crystal_String,
lambda:_?(NumericQ[#]||QuantityQ[#]&):-1,
condition___Condition,
OptionsPattern[]]:=Block[{
\[Lambda],limit,H,\[Theta],checkIfEmpty,custom,n,list,
G,Ginv,CrystalDot,res,
keep,options,angleThreshold,
(* Progress *)
progress,total
},

(* Dynamical progress *)
	progress={0,"Initialisation"};
	total=11;
	If[OptionValue["ShowProgress"],
	PrintTemporary[Row[{
	ProgressIndicator@Dynamic[progress[[1]]/total],
	Spacer[10],
	Dynamic@progress[[2]]
	}]]
	];

(* Checking input *)
	progress={0,"Checking input"};
	\[Lambda]=InputCheck["ProcessWavelength",crystal,lambda];

	limit=OptionValue["Limit"];
	If[!(Positive[limit]&&IntegerQ[limit]),
	Message[ReflectionList::limit];Abort[]];

	(* Useful variables *)
	progress={2,"Defining variables"};
	G=GetCrystalMetric[crystal];
	Ginv=Inverse@G;
	H=N@Chop@Ginv;
	\[Theta][hkl_]:=N[ArcSin[\[Lambda]*Sqrt[hkl.H.hkl]/2]];

		(* Empty list check *)
		checkIfEmpty:=If[list=={},
		Message[ReflectionList::empty];
		Goto["ReturnEmpty"]];

	(* Option: Use custom reflections *)
	progress={3,"Custom input"};
	If[ListQ@OptionValue["CustomReflections"],
	custom=OptionValue["CustomReflections"];

	(* Check custom input *)
	If[Flatten@custom=={},
	Message[InputCheck::hkl];
	Return[{}]
	];

	Check[InputCheck[custom,"Integer","WrapSingle"],
	Abort[]];	

	list=custom;
	Goto["ListDone"]
	];

(** Generating a reflection list **)
	progress={4,"Generating a reflection list"};
	(* Dot product in reciprocal space *)
	CrystalDot[u_,v_]:=Return[
	Sum[
	Sum[
	Ginv[[i,j]]*u[[i]]*v[[j]],
	{j,3}],
	{i,3}]
	];

	(* Coarse decision on which 'hkl' values to generate *)
	progress={5,"Deciding which 'hkl' values to generate"};
	options=#->OptionValue[#]&/@Keys@Options@ReflectionList;
	n=1;
	While[Im@\[Theta][{n,n,n}]==0,n++];
	n=Min[n,limit];
	list=ReflectionList[n,condition,options];
	checkIfEmpty;

	(* Filter away reflections with complex Bragg angle *)
	progress={6,"Checking Bragg angles"} ;
	list=Select[list,
	Norm[#]<=1&&Head[#]=!=Complex&[\[Theta][#]]&];

	(* Optional: Truncate at chosen angle threshold *)
	angleThreshold=OptionValue["AngleThreshold"];
	If[0<=angleThreshold<\[Pi]/2,
	list=Select[list,(\[Theta][#]<=angleThreshold)&]
	];

	(* Resolution filtering *)
	progress={7,"Resolution filtering"};
	res=\[Lambda]/2;
	list=Select[list,Sqrt[CrystalDot[#,#]]<1/(1.01*res)&];
	checkIfEmpty;

	(* Filter away absent reflections *)
	Label["ListDone"];
	progress={8,"Filtering away absent reflections"};
	list=Pick[
	list,
	SystematicAbsentQ[crystal,list],
	False];
	checkIfEmpty;


(** Optional: Merge symmetry-equivalent reflections **)
	progress={9,"Merging symmetry equivalent reflections"};
	If[OptionValue["SplitEquivalent"],
	(* Split *)
	Null,
	(* Merge *)
	list=MergeSymmetryEquivalentReflections[
	crystal,list,
	"ToStandardSetting"->OptionValue["ToStandardSetting"]]
	];

(** Optional: Limit number of reflections **)
	progress={10,"Limiting the number of reflections"};
	keep=OptionValue["Keep"];
		(* Check *)
		If[keep==All,keep=Length@list];
		If[!(IntegerQ[keep]&&Positive[keep]),
		Message[ReflectionList::keep]];


	progress={11,"Reflection list done"};
	list=Take[list,Min[keep,Length@list]];

(*---* End *---*)
Return@list;

Label["ReturnEmpty"];
Return[{}]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
Options@RelatedFunctionsGraph={
"Limit"->10,
"DirectOnly"->False,
"ShowDependent"->False
};

SetAttributes[RelatedFunctionsGraph,HoldFirst];

SyntaxInformation@RelatedFunctionsGraph={
"ArgumentsPattern"->{_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
RelatedFunctionsGraph[function_,OptionsPattern[]]:=Block[{
f,allf,deffile,import,
finddepf,data,d,main,g,done,new,x,X,c,part},

(* Loading data *)
f=ToString@HoldForm@function;
allf=First/@DeleteCases[Flatten@First@$MaXrdFunctions,""];
deffile=FileNameJoin[{$MaXrdPath,"Core","Definitions.m"}];
import=StringJoin@Check[Import[deffile,"Text"],Abort[]];

	(* Optional: Show all functions depedent on 'f' *)
	If[OptionValue["ShowDependent"],
	data=StringCases[import,Shortest[
	c:(allf)~~{"[",":="}~~d__~~"End[];"]:>{c,d}];
	x=Reap@Do[
	c=data[[i,1]];
	d=Sort@DeleteDuplicates@StringCases[data[[i,2]],allf];
	d=DeleteCases[d,c];
	Sow[c->d],
	{i,Length@data}];
	x=x[[2,1]];
	part=First/@Position[x,f];
	g=First/@x[[part]];
	g=DeleteCases[g,f];
		(* If none, return empty list *)
		If[g=={},Return[{}]];
	g=f->#&/@g;
		Goto["GraphDataGenerated"]
	];

(* Function for finding related functions *)
finddepf[F_]:=(
data=StringCases[
import,Shortest[
"Begin[\"`Private`\"];"~~Whitespace~~
"(* ::Input::Initialization:: *)\n"~~
F~~{"[",":="}~~d__~~"End[];"]
:>d];

	(* Check if anything is found *)
	If[data=={},Return[{}]];

d=DeleteDuplicates@Flatten@Sort@StringCases[First@data,{
"\""~~allf~~"\"",
allf~~"::",
allf}];
d=DeleteCases[d,F];
d=DeleteCases[d,x_/;StringContainsQ[x,{"\"","::"}]]
);

(* Seed *)
main=finddepf@f;
	If[main=={},Return[{}]];

(* Loop *)
g={};
done={};
new={f};
While[new!={},
x=First@new;
X=finddepf@x;
g=DeleteDuplicates@Join[g,#->x&/@X];

AppendTo[done,x];
new=Join[new,X];
new=Complement[new,done]
];

Label["GraphDataGenerated"];
g=DeleteDuplicatesBy[g,Sort];

	(* Optional: Only directly connected *)
	If[OptionValue["DirectOnly"],
	g=DeleteCases[g,(a_->b_)/;(a!=f&&b!=f)]];
	
	(* Optional: Limiting graph *)
	g=Take[g,UpTo@OptionValue["Limit"]];

(* Plot *)
(* A. Older than version 12 *)
If[$VersionNumber<12.,
GraphPlot[g,
DirectedEdges->True,
ImageSize->Large,
MultiedgeStyle->False,
DirectedEdges->True,
SelfLoopStyle->None,
VertexLabeling->True,
VertexRenderingFunction->({
White,EdgeForm[],Rectangle[
#-{0.02*StringLength@#2,0.05},#+{0.02*StringLength@#2,0.05}],Black,
Text[
Style[Hyperlink[#2, "paclet:MaXrd/ref/"<>#2],
11,"Program"],
#1]
}&),
Method->"RadialDrawing"],

(* B. Version 12 and above *)
GraphPlot[g,
DirectedEdges->True,
VertexLabels->None,
ImageSize->Large,
MultiedgeStyle->False,
DirectedEdges->True,
EdgeShapeFunction->({Arrowheads[{{Automatic,0.5}}],Arrow[#1]}&),
SelfLoopStyle->None,
VertexShapeFunction->({
White,EdgeForm[],Rectangle[
#-{0.02*StringLength@#2,0.05},#+{0.02*StringLength@#2,0.05}],Black,
Text[
Style[Hyperlink[#2, "paclet:MaXrd/ref/"<>#2],
11,"Program"],
#1]
}&),
Method->"RadialDrawing"]
]
];


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
ResetCrystalData::DemoDataNotFound="Demo data not found.";

SyntaxInformation@ResetCrystalData={
"ArgumentsPattern"->{}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
ResetCrystalData[]:=Block[{demoFile,newDataFile},
demoFile=FileNameJoin[
{$MaXrdPath,"Core","Data","CrystalDataDemo.m"}];
If[!FileExistsQ@demoFile,
Message[ResetCrystalData::DemoDataNotFound];Return[]];
newDataFile=FileNameJoin[{$MaXrdPath,"UserData","CrystalData.m"}];
CopyFile[demoFile,newDataFile,OverwriteTarget->True];
$CrystalData=Import@newDataFile;
Keys@$CrystalData
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
SimulateDiffractionPattern::InvalidStructureInput="Structural input must be a crystal label or the path to one or more structure files.";
SimulateDiffractionPattern::InvalidInputForDIFFUSE="If not using a crystal label with DIFFUSE, input should be two structure files.";
SimulateDiffractionPattern::InvalidReciprocalPlane="Invalid reciprocal plane input.";
SimulateDiffractionPattern::InvalidReciprocalSpaceLimit="Invalid setting for \"IndicesLimit\".";
SimulateDiffractionPattern::ZeroIntensity="No intensity found in data.";
SimulateDiffractionPattern::MissingOutputData="Unable to import the expected output data.";
SimulateDiffractionPattern::MissingProgram="`1` does not appear to be installed.";
SimulateDiffractionPattern::InvalidPrint="Invalid print setting.";
SimulateDiffractionPattern::InvalidFormat="Structure file seems to be invalid.";
SimulateDiffractionPattern::UnsupportedProgram="The program \[LeftGuillemet]`1`\[RightGuillemet] is not supported.";
SimulateDiffractionPattern::InvalidSubtractionMode="Invalid scattering subtraction mode.";
SimulateDiffractionPattern::InvalidScalingFactor="The scaling factor must be a number.";

Options@SimulateDiffractionPattern=SortBy[Normal@Union[
Association@Options@ArrayPlot,
Association@Options@ListDensityPlot,<|
"BraggScatteringSubtractionMode"->None,
"IndicesLimit"->5.5,
"LowerCutoff"->0,
"PrintOutput"->"ErrorsOnly",
"ProgramPaths"-><|
"MacOSX"-><|
"DISCUS"->"/usr/local/bin",
"DIFFUSE"->"/usr/local/bin"
|>,
"Windows"-><|
"DISCUS"->"C:\\Program Files (x86)\\Discus\\bin\\discus.exe",
"DIFFUSE"->""
|>
|>,
"ScalingFactor"->1,
"UseRawInput"->False,
(* ArrayPlot *)
ColorFunction->"Warm",
Frame->False,
FrameTicks->All,
ImageSize->Large,
PlotLegends->None,
ScalingFunctions->"Log"
|>],ToString[#[[1]]]&];

SyntaxInformation@SimulateDiffractionPattern={
"ArgumentsPattern"->{_,_,OptionsPattern[
{SimulateDiffractionPattern,ArrayPlot,ListDensityPlot}]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
SimulateDiffractionPattern[usingProgram_String,structureInput_,ImagePlane_,OptionsPattern[]]:=Block[{
imgPlane=ImagePlane,programPaths,searchExpression,options,inputs
},

(*---* Common driver routine *---*)
(* Common checks *)
If[!MemberQ[{"DISCUS","DIFFUSE"},usingProgram],
Message[SimulateDiffractionPattern::UnsupportedProgram,usingProgram];Abort[]];

If[!MemberQ[{"ErrorsOnly",All},OptionValue["PrintOutput"]],
Message[SimulateDiffractionPattern::InvalidPrint];Abort[]];

If[!MemberQ[
{None,"Biso","ExactAverage","SmallAverage"},
OptionValue["BraggScatteringSubtractionMode"]],
Message[SimulateDiffractionPattern::InvalidSubtractionMode];Abort[]];

If[!NumericQ@OptionValue["ScalingFactor"],
Message[SimulateDiffractionPattern::InvalidScalingFactor];Abort[]];

Which[
usingProgram==="DISCUS",
If[!StringQ@structureInput,
Message[SimulateDiffractionPattern::InvalidStructureInput];
Abort[]],

usingProgram==="DIFFUSE",
If[StringQ@structureInput,InputCheck["CrystalQ",structureInput],
If[Length@structureInput!=2||AnyTrue[structureInput,!FileExistsQ[#]&],
Message[SimulateDiffractionPattern::InvalidInputForDIFFUSE];
Abort[]]
]
];

If[StringQ@imgPlane,imgPlane=MillerNotationToList@imgPlane];
If[MatchQ[Sort@imgPlane,{_Integer,#,#}]
&["h"|"k"|"l"]\[Nand]DuplicateFreeQ@imgPlane,
Message[SimulateDiffractionPattern::InvalidReciprocalPlane];Abort[]];

If[!DirectoryQ@#,CreateDirectory@#]&[
FileNameJoin[{$TemporaryDirectory,"MaXrd"}]];

programPaths=OptionValue["ProgramPaths"][$OperatingSystem][usingProgram];
If[programPaths===""||DirectoryQ@programPaths,
searchExpression=Which[
usingProgram==="DISCUS",{"discus"},
usingProgram==="DIFFUSE",{"dzmc","bin2gray"}
];
If[$OperatingSystem==="Windows",searchExpression=#<>".exe"&/@searchExpression];
searchExpression="(?i)"<>StringRiffle[searchExpression,"|"];
programPaths=Sort@FileNames[
RegularExpression@searchExpression,programPaths,IgnoreCase->True];
If[programPaths==={}||!AllTrue[programPaths,FileExistsQ],
Message[SimulateDiffractionPattern::MissingProgram,usingProgram];
Abort[]];
];
If[usingProgram==="DISCUS"&&ListQ@programPaths,
programPaths=First@programPaths];


(* Switch flow *)
options=Thread[#->OptionValue[#],String]&/@Keys@Options@SimulateDiffractionPattern;
inputs={programPaths,structureInput,imgPlane,Association@options};
Which[
usingProgram==="DISCUS",SDP$DISCUS@@inputs,
usingProgram==="DIFFUSE",SDP$DIFFUSE@@inputs
]
]


(* ::Input::Initialization:: *)
SDP$DISCUS[programPath_String,structureInput_,ImagePlane_,givenOptions_]:=Block[{
structureFile=structureInput,options=givenOptions,workDir,
i,stream,line,streamData,allCoords,streamLength,
latticeParameters,crystalM,
structureSize,sizeX,
hklMax,abscissaIndex,ordinateIndex,
x,imageOrientation,commands,feedback="",
cutOffValue,data,dataLength,
xDataSize,yDataSize,xMin,xMax,yMin,yMax,xStep,yStep,numbers,plotData,
scaleMax,intensities,maxIntensity,useRawInputQ,imageBasis
},

(* Handle both crystal label and structure file input *)
If[KeyExistsQ[$CrystalData,structureInput],
structureFile=ExportCrystalData["DISCUS",structureFile,FileNameJoin[
{$TemporaryDirectory,"MaXrd","TemporaryStructureFile.stru"}],
"IncludeStructureSizeInfo"->False
]];
If[!FileExistsQ@structureFile,
Message[SimulateDiffractionPattern::InvalidStructureInput];Abort[]];

(* Determining structure size *)
{i,stream}={1,OpenRead@structureFile};
While[True,
line=Read[stream,String];
If[StringTake[line,;;4]==="cell",
latticeParameters=line];
If[StringTake[line,;;5]==="atoms",
Break[]];
If[i>10,Message[SimulateDiffractionPattern::InvalidFormat];Abort[]];
i++;
];
streamData=ReadList[stream,String];
Close@stream;

allCoords=ToExpression@Part[StringSplit@streamData,All,2;;4];
streamLength=Length@allCoords;
structureSize=MinMax@allCoords;
sizeX=Round@structureSize[[2]];
If[structureSize[[1]]<-0.5,sizeX*=2];
sizeX=ToString@sizeX;
latticeParameters=ToExpression@StringCases[latticeParameters,{DigitCharacter,"."}..];
crystalM=GetCrystalMetric[latticeParameters,
"Space"->"Reciprocal","ToCartesian"->True];

(* Preparing input for Fourier transform *)
workDir=DirectoryName@structureFile;
hklMax=options["IndicesLimit"];
If[NumericQ@hklMax\[Nand]Positive@hklMax,
Message[SimulateDiffractionPattern::InvalidReciprocalSpaceLimit];Abort[]];

imageOrientation=InputCheck["GetReciprocalImageOrientation",
latticeParameters,ImagePlane,hklMax,False];

{abscissaIndex,ordinateIndex}={#1,#2}&@@Flatten[
Position[ImagePlane,#]&/@{"h","k","l",_Integer}];

commands="
cd "<>workDir<>"
################################################
# COMBINED BUILD MACRO FOR `SimulateDiffractionPattern`
################################################
   reset
####### Load/build crystal #####################
variable int, sizeX
sizeX = "<>sizeX<>"
#
read
  stru "<>FileNameTake@structureFile<>"
#
chem
  elem
exit
####### Fourier transform ######################
variable real, hklMax
variable int,  fourierWidth
variable int,  fourierPoints
#
hklMax = "<>ToString@N@hklMax<>"
fourierWidth = 2 * hklMax
fourierPoints = fourierWidth * sizeX + 1
#
fourier
  xray
  wvle moa1
#
  ll   "<>imageOrientation[[1]]<>"
  lr   "<>imageOrientation[[2]]<>"
  ul   "<>imageOrientation[[3]]<>"
  na   fourierPoints
  no   fourierPoints
  abs  "<>(abscissaIndex/.{1->"h",2->"k",3->"l"})<>"
  ord  "<>(ordinateIndex/.{1->"h",2->"k",3->"l"})<>"
#
  show
  run
exit
#
#
#---# Fourier data output #---#
output
  value intensity
  format standard
  outfile fourier_data.dat
  run
exit
################################################
exit
";

(* Run DISCUS *)
feedback=RunProcess[programPath,All,commands];
SDP$EvaluateFeedbackPrint[commands,feedback,options["PrintOutput"]];

(*-----* Plot preparations *-----*)
(* Importing (x,y,intensity) data from file *)
data=Check[
Import[FileNameJoin[{workDir,"fourier_data.dat"}],"Table"],
Message[SimulateDiffractionPattern::MissingOutputData];Abort[]];
dataLength=Length@data;

i=1;
While[data[[i,1]]==="#"&&i<=dataLength,i++];
Check[{xDataSize,yDataSize}=data[[i]],Abort[]];
{xMin,xMax,yMin,yMax}=data[[i+1]];
xStep=(xMax-xMin)/(xDataSize-1);
yStep=(yMax-yMin)/(yDataSize-1);

numbers=Flatten[data[[i+2;;]]];
numbers=Partition[numbers,xDataSize];

plotData=Table[
{
xMin+(x-1)*xStep,
yMin+(y-1)*yStep,
numbers[[y,x]](* Instead of transposing *)
},{y,yDataSize},{x,xDataSize}];
plotData=Flatten[plotData,1];

(* Scaling intensities *)
scaleMax=100.;
intensities=plotData[[All,3]];
maxIntensity=Max@intensities;
If[maxIntensity==0,Message[SimulateDiffractionPattern::ZeroIntensity];Abort[]];
intensities*=scaleMax/maxIntensity;
intensities=#*options["ScalingFactor"]&/@intensities;
intensities=intensities/.x_/;x>scaleMax->scaleMax;(* 16 bit max *)
plotData[[All,3]]=intensities;

(* Data treatment and preparation *)
useRawInputQ=TrueQ@options["UseRawInput"];
If[useRawInputQ,
(* a. Use data "as is" *)
plotData=Partition[plotData[[All,3]],Length@numbers],

(* b. Rescale intensity and use appropriate basis for image *)
cutOffValue=Power[10.,-3];
plotData=plotData/.{x_,y_,i_}/;i<cutOffValue:>{x,y,cutOffValue};
imageBasis=Normalize/@crystalM[[
{abscissaIndex,ordinateIndex},
{abscissaIndex,ordinateIndex}]];
plotData[[All,{1,2}]]=Map[imageBasis.#&,plotData[[All,{1,2}]]]
];

(* Prepare and deliver plot *)
If[useRawInputQ,
AssociateTo[options,DataRange->{{xMin,xMax},{yMin,yMax}}],
AssociateTo[options,AspectRatio->Divide@@Total@imageBasis]
];

ListDensityPlot[plotData,Sequence@@FilterRules[Normal@options,Options@ListDensityPlot]]
]


(* ::Input::Initialization:: *)
SDP$DIFFUSE[programPaths_List,structureInput_,ImagePlane_List,givenOptions_]:=Block[{
structureFiles,workDir=FileNameJoin[{$TemporaryDirectory,"MaXrd"}],
options=givenOptions,inputFileDZMC,
subtractionMode,lowerCutoff,
inputFile1="diffuse_input1_crystal.txt",
inputFile2="diffuse_input2_setup.txt",
commands,feedback,outputFile,imageData
},

(* Handle both crystal label and structure file input *)
If[KeyExistsQ[$CrystalData,structureInput],
subtractionMode=options["BraggScatteringSubtractionMode"];
subtractionMode=subtractionMode/.{
None->"N","Biso"->"Y","ExactAverage"->"E","SmallAverage"->"e"};
structureFiles=ExportCrystalData["DIFFUSE",structureInput,
workDir,ImagePlane,options["IndicesLimit"],subtractionMode],

CopyFile[#1,#2,OverwriteTarget->True]&@@@Transpose[{structureInput,
FileNameJoin[{workDir,#}]&/@{inputFile1,inputFile2}
}]
];

If[AnyTrue[structureFiles,!FileExistsQ[#]&],
Message[SimulateDiffractionPattern::InvalidStructureInput];Abort[]];

Quiet@DeleteFile@FileNames["output*",workDir];
inputFileDZMC=FileNameJoin[{workDir,"DZMC_inputs.txt"}];
Put[OutputForm@inputFile2,OutputForm["output.bin"],inputFileDZMC];

commands=StringTemplate["
cd `workDir`
\"`prog1`\" `inp1` < DZMC_inputs.txt
\"`prog2`\" --quiet=true output.bin
"]@<|"workDir"->workDir,"inp1"->inputFile1,"prog1"->programPaths[[2]],"prog2"->programPaths[[1]]|>;

(* Run DIFFUSE and then bin2gray *)
feedback=RunProcess[$SystemShell,All,commands];
SDP$EvaluateFeedbackPrint[commands,feedback,options["PrintOutput"]];

outputFile=FileNameJoin[{workDir,"output.pgm"}];
If[!FileExistsQ@outputFile,
Message[SimulateDiffractionPattern::MissingOutputData];Abort[]];

lowerCutoff=options["LowerCutoff"];
If[TrueQ@options["UseRawInput"],
Import@outputFile,

imageData=N@Import[outputFile,"Data"];
imageData=#*options["ScalingFactor"]&/@imageData;
imageData=imageData/.x_/;x>65535->65535;(* 16 bit max *)
imageData=imageData/.x_/;x<lowerCutoff->0.;
ArrayPlot[imageData,Sequence@@FilterRules[Normal@options,Options@ArrayPlot]]
]
]


(* ::Input::Initialization:: *)
SDP$EvaluateFeedbackPrint[commands_String,Feedback_Association,optionSetting_]:=Block[{feedback=Feedback,stderr},
If[optionSetting===All,
Print@Prepend[feedback,"Input"->commands],

(* Removing irrelevant errors *)
stderr=StringTrim@StringDelete[feedback["StandardError"],{
"Remaining memory"~~__~~WhitespaceCharacter,
WhitespaceCharacter~~"More segments"~~__,
"'\\\\"~~__~~"UNC paths"~~__~~"Defaulting to Windows directory."
}];
If[stderr=!=""||feedback["ErrorCode"]==2,
Print@stderr]
]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
StructureFactor::invalidsetting="Invalid setting of the `1` option.";
StructureFactor::extinct="The reflection `1` is systematically absent for space group `2`.";
StructureFactor::mismatch="Element mismatch detected.";

Options@StructureFactor={
"AbsoluteValue"->True,
"Units"->True,
"IgnoreSystematicAbsence"->False,
"Threshold"->Power[10.,-6],
(* Options from 'GetAtomicScatteringFactors' *)
"DispersionCorrections"->OptionValue[
GetAtomicScatteringFactors,"DispersionCorrections"],
"f0Source"->OptionValue[
GetAtomicScatteringFactors,"f0Source"],
"f1f2Source"->OptionValue[
GetAtomicScatteringFactors,"f1f2Source"],
(* Options from 'GetElements' *)
"IgnoreIonCharge"->OptionValue[
GetElements,"IgnoreIonCharge"]
};

SyntaxInformation@StructureFactor={
"ArgumentsPattern"->{_,_,_.,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
StructureFactor[
crystal_String,
hklInput_List,
lambda:_?(NumericQ[#]||QuantityQ[#]&):-1,
OptionsPattern[]]:=Block[{
data,abs,ignoreExtinct,units,\[Delta],hkl,L0,L,\[Lambda],
zerotype,absence,l,hklPos,j,
sgHM,H,d,sl,fOptions,f,
atomdata,r,operators,type,
disp,U,R,T,lattice,cvecs,occ,elements,
siteSymMxyz,siteSymO,
sf,SF,F,\[Phi]},

(*---* Checking input *---*)
	(* Crystal and wavelength/energy *)
	InputCheck["CrystalQ",crystal];
	\[Lambda]=InputCheck["ProcessWavelength",crystal,lambda];
	data=$CrystalData[crystal];

	(* Option settings *);
	abs=OptionValue["AbsoluteValue"];
	ignoreExtinct=OptionValue["IgnoreSystematicAbsence"];
	units=OptionValue["Units"];
	\[Delta]=OptionValue["Threshold"];
		If[!NumericQ@\[Delta],
		Message[StructureFactor::"invalidsetting","\"Threshold\""];
		Abort[]];

	hkl=InputCheck[hklInput,"Integer","WrapSingle"];
	L0=Length@hkl;

(*---* Systematically absent reflections *---*)
	If[!ignoreExtinct,
	(* a. Check for systematic absence *)
	zerotype=If[abs,
	If[units,
	{0,Quantity[0,"Degrees"]},
	{0,0}],
	0];

	absence=SystematicAbsentQ[crystal,hkl];
	l=Length@absence;
	hkl=Pick[hkl,absence,False];
	hklPos=Position[absence,False];

	(* Return if only extinct reflections *)
	If[AllTrue[Flatten@absence,TrueQ],
	If[l===1,
	Return@zerotype,
	Return@ConstantArray[zerotype,l]
	]],

	(* b. Keep list as is *)
	Null
	];

(*---* Useful variables *---*)
	(* Miscellaneous *)
	L=Length@hkl;
	sgHM=GetSymmetryData[data[["SpaceGroup"]],"HermannMauguinFull"];
	H=Chop@N@Inverse@GetCrystalMetric[crystal];
	sl=N[Sqrt[#.H.#]/2]&/@hkl;

	(* Obtaining atomic scattering factors *)
	fOptions=Keys@Options@GetAtomicScatteringFactors;
	fOptions=DeleteCases[fOptions,"SeparateCorrections"];
	fOptions=#->OptionValue[#]&/@fOptions;

	f=GetAtomicScatteringFactors[crystal,hkl,\[Lambda],Sequence@@fOptions];
	(* Wrap 'f' if single reflection input *)
	If[AssociationQ@f,f={f}];

	atomdata=data["AtomData"];
	r=atomdata[[All,"FractionalCoordinates"]];

	(* Variables related to displacement parameters *)
	d=Chop@N@Sqrt@DiagonalMatrix@Diagonal@H;
	disp=Lookup[atomdata,"DisplacementParameters",0.];
	
	(* Atomic displacement parameters preparation *)
	U={};
	Do[
	If[Length@disp[[i]]==6,
	AppendTo[U,
	Partition[Part[disp[[i]],{1,4,5,4,2,6,5,6,3}],3]],
	AppendTo[U,disp[[i]]]
	],
	{i,Length@disp}];

	(* Symmetry operators *)
	operators=GetSymmetryOperations[crystal];
	{R,T}=Transpose@operators;
	
	elements=atomdata[[All,"Element"]];
		(* Optional: Remove charge/ions *)
		If[OptionValue["IgnoreIonCharge"],
		elements=StringDelete[elements,
		{DigitCharacter,"+","-"}]
		];
		(* Check if elements don't match 'f' *)
		If[Sort@Keys@First@f=!=Sort@DeleteDuplicates@elements,
		Message[StructureFactor::mismatch];Abort[]];

	type=Lookup[atomdata,"Type","Uiso"];

	(* Centring *)
	lattice=StringTake[sgHM,1];
	cvecs=Length@InputCheck["GetCentringVectors",lattice];

	(* Site symmetry multiplicity and order *)
		(* Site multiplicity of general position *)
		siteSymMxyz=cvecs*Length@operators;
		(* Site symmetry order: siteSymO = siteSymMxyz / siteSymM *)

	Which[
	(* a. Extract 'SiteSymmetryOrder' from $CrystalData *)
	KeyExistsQ[First@atomdata,"SiteSymmetryOrder"],
	siteSymO=atomdata[[All,"SiteSymmetryOrder"]],

	(* b. Calculate order from 'SiteSymmetryMultiplicity' *)
	KeyExistsQ[First@atomdata,"SiteSymmetryMultiplicity"],
	siteSymO=siteSymMxyz/
	atomdata[[All,"SiteSymmetryMultiplicity"]],

	(* c. Calculate site symmetry order *)
	True,
	siteSymO=siteSymMxyz/Table[Length@
	SymmetryEquivalentPositions[
	crystal,r[[a]]],{a,Length@r}]
	];

	(* Occupation factor *)
	Label["OccupationFactor"];
	occ=Lookup[atomdata,"OccupationFactor",1.];


(*---* Structure factor calculation *---*)
(* Magnitude *)
	sf=Table[(* Table for each reflection *)
		Sum[1
	*occ[[j]](* Occupation factor *)
	*cvecs(* Centring vectors *)
	*1.0/siteSymO[[j]] (* Symmetry reduction *)
	*Part[f[[h]],elements[[j]]](* Atomic form factor *)
	*Sum[
	Which[(* Atomic displacement *)
	type[[j]]=="Uani",
		Exp[-2 Pi^2*hkl[[h]].d.R[[s]]
		.U[[j]]
		.Transpose[R[[s]]]
		.d.hkl[[h]]],
	type[[j]]=="Uiso",
		Exp[-8 Pi^2*disp[[j]]*(sl[[h]])^2],
	type[[j]]=="Bani",
		Exp[-1/4*hkl[[h]].d.R[[s]]
		.U[[j]]
		.Transpose[R[[s]]]
		.d.hkl[[h]]],
	type[[j]]=="Biso",(* Temperature factor *)
		Exp[-disp[[j]]*(sl[[h]])^2],
	True,Message[$CrystalData::type,type];Abort[]
	]*Exp[2Pi*I(hkl[[h]].R[[s]].r[[j]]+hkl[[h]].T[[s]])],
	{s,Length@R}],
		{j,Length@atomdata}],
			{h,L}];

(* Phase *)
	If[!abs,SF=sf;Goto["ComplexNumber"]];

	SF=Reap[Do[
	F=sf[[i]];
	If[
	(* a. Check threshold *)
	Abs@Re@F<\[Delta],Sow[{0,0}],
	(* b. Calculate *)
	Which[
	Re[F]==0.&&Im[F]>=0,
		F=0;\[Phi]=90,
	Re[F]==0.&&Im[F]<0,
		sf[[i]]=0;\[Phi]=-90,
	Re[F]>=0.&&Im[F]==0,
		\[Phi]=90,
	Re[F]<0.&&Im[F]==0,
		\[Phi]=180,
	Re[F]>0.&&Im[F]>0,	
		\[Phi]=N[ArcTan[Abs[Im[F]/Re[F]]]/Degree],
	Re[F]>0.&&Im[F]<0,
		\[Phi]=-N[ArcTan[Abs[Im[F]/Re[F]]]/Degree],
	Re[F]<0.&&Im[F]>0,
		\[Phi]=N[(Pi-ArcTan[Abs[Im[F]/Re[F]]])/Degree],
	Re[F]<0.&&Im[F]<0,
		\[Phi]=N[(ArcTan[Abs[Im[F]/Re[F]]]-Pi)/Degree]
	];
	Sow[{F,\[Phi]}]],
	{i,L}]][[2,1]];


(*---* Preparing output *---*)
	(* Processing units *)
	If[units,SF=MapAt[Quantity[#,"Degrees"]&,SF,{All,2}]];

	(* Organising structure factors *)
		(* a. Return absolute value and phase *)
		SF=MapAt[Abs,SF,{All,1}];

		(* b. Return complex number *)
		Label["ComplexNumber"];

	(* Putting back extinct reflections *)
	If[L0>1&&!ignoreExtinct,SF=ReplacePart[
	ConstantArray[zerotype,l],
	Thread[hklPos->SF]]];

If[L0===1,First@SF,SF]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
StructureFactorTable::invcustom="Invalid setting for the \[LeftGuillemet]CustomReflections\[RightGuillemet] option.";
StructureFactorTable::sort="Invalid \[LeftGuillemet]Sort\[RightGuillemet] option.";

Options[StructureFactorTable]={
(* DarwinWidth and ExtinctionLength options *)
"Polarisation"->"\[Pi]",
(* StructureFactor options *)
"Threshold"->Power[10.,-6],
(* ReflectionList options *)
"SplitEquivalent"->False,
"CustomReflections"->False,
"ReflectionListKeep"->50,
(* StructureFactorTable options *)
"Keep"->All,
"Sort"->-2,
"TitleStyle"->{FontFamily->"Baskerville",FontSize->15},
"SubtitleStyle"->{FontFamily->"Times New Roman",FontSize->13,Gray},
"NumberStyle"->FontFamily->"Courier",
Background->{{None},{None,{None,LightGray}}},
Dividers->{{None,{True},None},{None,None,True}}
};

SyntaxInformation@StructureFactorTable={
"ArgumentsPattern"->{_,_.,_.,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
StructureFactorTable[
crystal_String,
lambda:_?(NumericQ[#]||QuantityQ[#]&):-1,
condition___Condition,
OptionsPattern[]]:=
Block[{
\[Lambda],hkl,H,sl,bragg,
column1,column2,column3,column4,column5,column6,
V,
zeros,
R,\[Theta],Fh,FhBar,\[CapitalLambda]o,\[Delta]os,
temp,temp1,temp2,temp3,temp4,temp5,temp6,
sort,keep,polarisation,threshold},

(*---* Preparations *---*)
	(* Checking wavelength *)
	\[Lambda]=InputCheck["ProcessWavelength",crystal,lambda];

	(* Preparing list of reflections *)
	hkl=Check[ReflectionList[
	crystal,\[Lambda],condition,
	"SplitEquivalent"->OptionValue["SplitEquivalent"],
	"CustomReflections"->OptionValue["CustomReflections"],
	"Keep"->OptionValue["ReflectionListKeep"]
	],Abort[]];

	(* Useful variables *)
	H=Chop@N@Inverse@GetCrystalMetric[crystal];
	(* Using inner product and Bragg's law *)
	sl[h_]:=N[Sqrt[h.H.h]/2];
	bragg[h_]:=N[ArcSin[sl[h]*QuantityMagnitude@\[Lambda]]];

	(* Miscellaneous options *)
	polarisation=OptionValue["Polarisation"];
	threshold=OptionValue["Threshold"];

(*---* Calculations *---*)
	(* Miller indices *)
	column1=MillerNotationToString/@hkl;

	(* Structure factors *)
	temp=Check[StructureFactor[crystal,hkl,\[Lambda],"Threshold"->threshold],Abort[]];
	column2=Chop[temp[[All,1]],OptionValue["Threshold"]];
		(* Structure factor zero positions *)
		zeros=Position[column2,0];

	(* Phase *)
	column3=QuantityMagnitude@temp[[All,2]];
	column3=ReplacePart[column3,zeros->"\[Dash]"];

	(* Bragg angles *)
	column4=(bragg/@hkl)/Degree;
	
	(* Unit cell volume *)
	V=Sqrt@Det@GetCrystalMetric[crystal];

(* Extinction length (Pendell\[ODoubleDot]sung distance) *)
	R=QuantityMagnitude@UnitConvert[
	Quantity["ClassicalElectronRadius"],"Angstroms"];

	\[Theta]=column4*Degree;
	Fh=column2;
	FhBar=(StructureFactor[crystal,-hkl,\[Lambda]])[[All,1]];
	
	\[CapitalLambda]o=Quiet@ExtinctionLength[crystal,\[Lambda],hkl,
"Units"->False,"Polarisation"->polarisation];
	\[CapitalLambda]o=ReplacePart[\[CapitalLambda]o,zeros->"\[Dash]"];
	
(* Darwin width *)
	\[Delta]os=Quiet@DarwinWidth[crystal,\[Lambda],hkl,
"Units"->False,"Polarisation"->polarisation];
	\[Delta]os=ReplacePart[\[Delta]os,zeros->"\[Dash]"];
	
	{column5,column6}={\[CapitalLambda]o,\[Delta]os};

(*---* Preparing output *---*)
	temp1={column1,column2,column3,column4,column5,column6};

	(* Optional: Sorting option *)
		sort=OptionValue["Sort"];
		
		(* Check sort setting *)
		If[!MemberQ[Range[6],Abs@sort],
		Message[StructureFactorTable::sort];
		Abort[]];
		
		(* Sorting by a specific column *)			
		temp2=Sort[Transpose@temp1,
		#1[[Abs@sort]]<#2[[Abs@sort]]&];
			(* Reversing if negative *)
			If[sort<0,temp2=Reverse@temp2];

	(* Optional: Truncate the list of reflections *)
		keep=OptionValue["Keep"];
		If[keep==All,keep=Length@hkl];
		If[!IntegerQ[keep]||!Positive[keep],
		Message[StructureFactorTable::invkeep];Abort[]];
		
		temp2=Take[temp2,Min[keep,Length@hkl]];

(* Rounding off numbers *)
	temp3=NumberForm[#,{5,3},
	DigitBlock->3,NumberSeparator->","]&/@(Flatten@temp2);
	temp4=Partition[temp3,Last@Dimensions@temp2];


(*---* Table construction *---*)
	temp5=PrependTo[temp4,
	{"(hkl)","|\!\(\*SubscriptBox[\(F\), \(hkl\)]\)|","\!\(\*SubscriptBox[\(\[Phi]\), \(hkl\)]\) [\[Degree]]",
	"\!\(\*SubscriptBox[\(\[Theta]\), \(B\)]\) [\[Degree]]","\!\(\*SubscriptBox[\(\[CapitalLambda]\), \(0\)]\) [\[Micro]m]","2\!\(\*SubscriptBox[\(\[Delta]\), \(os\)]\) [\[Micro]rad]"}];
	temp6=PrependTo[temp5,
	{Null,"Structure factor","Phase",
	"Bragg angle","Extinction length",
	"Darwin width"}];

	Grid[temp6,
	Dividers->OptionValue[Dividers],
	Background->OptionValue[Background],
	Alignment->{Center,Center},
	Spacings->{Automatic,2->0.5},
	ItemStyle->
	{Automatic,
	Automatic,
	{
	{{1,1},{1,First@Dimensions@temp1}}->OptionValue["TitleStyle"],
	{{2,2},{1,First@Dimensions@temp1}}->OptionValue["SubtitleStyle"],
	{{3,First@Dimensions@temp2+2},{1,First@Dimensions@temp1}}->OptionValue["NumberStyle"]
	}
	}]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
SymmetryEquivalentPositions::threshold="Tolerance specification must be a non-negative number.";

Options@SymmetryEquivalentPositions={
"UseCentring"->True,
"RationaliseThreshold"->0.001
};

SyntaxInformation@SymmetryEquivalentPositions={
"ArgumentsPattern"->{_,_.,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
SymmetryEquivalentPositions[input_String,xyzInput_List:{"x","y","z"},
OptionsPattern[]]:=
Block[{
group,\[Delta],useCentringQ,fracs,r,xyz,
s,R,T,equiv,rationalise,mod,generate,gather,centring,final,t,sym,c,add,pos,temp},

(*---* Input check *---*)
	group=InputCheck["GetPointSpaceGroupCrystal",input];
	xyz=InputCheck[xyzInput,"WrapSingle"];

	(* Options *)
	\[Delta]=OptionValue["RationaliseThreshold"];
	If[!(NumericQ[#]&&!Negative[#])&@\[Delta],
	Message@SymmetryEquivalentPositions::threshold];
	useCentringQ=OptionValue["UseCentring"];

(*---* Procedure for generating equivalent positions *---*)
(* Replace float zero to integer zero and rationalis *)
	fracs={1/12,1/8,1/6,1/4,1/3,3/8,5/12,1/2,7/12,5/8,2/3,3/4,5/6,7/8,11/12};
	rationalise[coord_]:=Reap[Do[Sow@If[
MemberQ[fracs,r=Rationalize[i/.{0.->0},\[Delta]]],r,i],
{i,coord}]][[2,1]];

(* Dealing with remainder of coordinates *)
	mod[X_]:=Which[
	NumericQ[X],Mod[X,1],
	Head[X]===Plus,
		t=Select[X,NumericQ];
		sym=X-t;
		Mod[t,1]+sym,
	True,X];

(* Preparing symmetry generation *)
	s=GetSymmetryOperations@group;
	{R,T}=Transpose@s;
	If[MatchQ[Dimensions@s,{_,3,3}],
	(* Point group procedure *)
	generate[coord_]:=Table[Transpose[s[[i]]].coord,{i,Length@s}],
	(* Space group procedure *)
	generate[coord_]:=Table[R[[i]].coord+T[[i]],{i,Length@s}]];
	
	(* Bring coordinates into one cell *)
	gather[list_]:=Map[mod,list,{2}];

(* Process centring *)
	c=StringTake[GetSymmetryData[group,"HermannMauguinFull"],1];
	add=InputCheck["GetCentringVectors",c];

	If[useCentringQ,
	(* a. Apply centring vectors *)
	centring[list_]:=(
	equiv=DeleteDuplicates@list;
	temp=equiv;
	equiv=Reap[Do[Sow[
	#+add[[i]]&/@temp],
	{i,Length@add}]][[2,1]];
	equiv=Flatten[equiv,1];
	Map[mod,equiv,{2}]
	),

	(* b. Do not use centring vectors *)
	centring[list_]:=(
	temp=Table[
	Map[mod,list[[i]]+add[[j]]],
	{j,Length@add},{i,Length@list}];

	pos=Table[Position[
list,#[[i]]],
{i,Length@list}]&/@temp;

	temp=Transpose@pos;
	temp=DeleteDuplicates@Table[
	Sort@Flatten[temp[[i]],1],
	{i,Length@temp}];

	temp=Flatten@temp[[All,1]];
	list[[Flatten@temp]]
	)
	];

	(* Final trimming *)
	final[list_]:=(
	temp=list/.x_Real:>Round[x,10.^(-6)];
	DeleteDuplicates@temp
	);


(*---* Execute procedure *---*)
xyz=final/@centring/@gather/@generate/@rationalise/@xyz;

If[MatchQ[xyzInput,{x_,y_,z_}/;!AnyTrue[{x,y,z},ListQ]],
Return@First@xyz,Return@xyz]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
SyntaxInformation@SymmetryEquivalentReflections={
"ArgumentsPattern"->{_,_.}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
SymmetryEquivalentReflections[group_String,hkl_List:{"h","k","l"}]:=Block[{s},
(* Input check *)
	If[!KeyExistsQ[$GroupSymbolRedirect,group]&&
	!KeyExistsQ[$CrystalData,group],
	Message[GetSymmetryOperations::missing,group];
	Abort[]];

	InputCheck[hkl,"1hkl","StringSymbol"];

(* Symmetry equivalent reflections *)
	If[MemberQ[Keys@$PointGroups,group],
	s=GetSymmetryOperations@group,
	s=GetSymmetryOperations@GetLaueClass@group];

	DeleteDuplicates@Table[
	Transpose[s[[i]]].hkl,{i,Length@s}]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
SyntaxInformation@SymmetryEquivalentReflectionsQ={
"ArgumentsPattern"->{_,_}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
SymmetryEquivalentReflectionsQ[group_String,hkl_List]:=
Block[{equiv},
(* Check input *)
Check[InputCheck[hkl,"Multiple"],Abort[]];

(* Listing all symmetry-equivalents of the first reflection *)
equiv=SymmetryEquivalentReflections[group,First@hkl];

(* Checking if all given reflections are symmetry equivalent *)
AllTrue[hkl,MemberQ[equiv,#]&]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
SynthesiseStructure::SizeError="Size discrepancy in domain input.";
SynthesiseStructure::DifferentBlockSizes="The blocks must have the same size.";
SynthesiseStructure::InvalidOutputSize="Output size must be a list of three positive integers.";
SynthesiseStructure::IncompatibleOutputSize="Output size must be compatible with block sizes.";
SynthesiseStructure::InvalidSelectionMethod="\[LeftGuillemet]`1`\[RightGuillemet] is not a valid selection method.";

Options@SynthesiseStructure={
"RotationAnchorReference"->"DomainCentroid",
"RotationAnchorShift"->{0,0,0},
"RotationAxes"->IdentityMatrix[3],
"RotationMap"-><||>,
"SelectionMethod"->"Sequential",
"UsePlacementBuffer"->False,
(* Options from 'EmbedStructure' *)
"OverlapPrecedence"->OptionValue[
EmbedStructure,"OverlapPrecedence"],
"OverlapRadius"->OptionValue[
EmbedStructure,"OverlapRadius"]
};

SyntaxInformation@SynthesiseStructure={
"ArgumentsPattern"->{_,_,_.,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
SynthesiseStructure[blocks_List,outputSize_List,outputName_String,OptionsPattern[]]:=Block[{
blockSizes,selectionMethod=OptionValue["SelectionMethod"],
insertionCoordinates,b,numberOfBlocks,blocksToUse,buildTasks
},
(* Checking input *)
Scan[InputCheck["CrystalQ",#]&,blocks];

If[!MemberQ[{"Random","Sequential"},selectionMethod],
Message[SynthesiseStructure::InvalidSelectionMethod,selectionMethod];
Abort[]];

(* Checking if all blocks have same size *)
blockSizes=($CrystalData[#,"Notes","StructureSize"]/._Missing->{1,1,1})&/@blocks;
If[!SameQ@@blockSizes,Message[SynthesiseStructure::DifferentBlockSizes];Abort[]];
blockSizes=blockSizes[[1]];

(* Checking if 'outputSize' is valid *)
If[(Length@outputSize!=3)||(AnyTrue[outputSize,Positive[#]\[Nand]IntegerQ[#]&]),
Message[SynthesiseStructure::InvalidOutputSize];Abort[]];

If[Total@MapThread[Mod,{outputSize,blockSizes}]=!=0,
Message[SynthesiseStructure::IncompatibleOutputSize];Abort[]];

(* Preparing placement of blocks *)
b=If[TrueQ@OptionValue["UsePlacementBuffer"],2,1];
insertionCoordinates=Flatten[Table[{i,j,k},
{i,0,b*outputSize[[1]]-1,b*blockSizes[[1]]},
{j,0,b*outputSize[[2]]-1,b*blockSizes[[2]]},
{k,0,b*outputSize[[3]]-1,b*blockSizes[[3]]}
],2];

numberOfBlocks=Times@@MapThread[Divide,{outputSize,blockSizes}];
blocksToUse=Flatten[ConstantArray[blocks,
\[LeftCeiling]numberOfBlocks/Length@blocks\[RightCeiling]]][[;;numberOfBlocks]];
If[selectionMethod==="Random",
blocksToUse=RandomSample@blocksToUse];

(* Assembling *)
buildTasks=Transpose[{blocksToUse,insertionCoordinates}];
AppendTo[$CrystalData,outputName->$CrystalData[buildTasks[[1,1]]]];
Scan[EmbedStructure[{#[[1]]},{#[[2]]},outputName,
"MatchHostSize"->False,"ShowProgress"->False,
"OverlapPrecedence"->OptionValue["OverlapPrecedence"],
"OverlapRadius"->OptionValue["OverlapRadius"]
]&,
buildTasks[[2;;]]];

$CrystalData[outputName,"Notes","StructureSize"]=Last@insertionCoordinates+blockSizes;

outputName
]


(* ::Input::Initialization:: *)
SynthesiseStructure[
{{A_Integer,B_Integer,C_Integer},domains_List},
outputName_String,
labelMap_:<||>,
OptionsPattern[]]:=Block[{
blocks,blockSizes,outputSize,b,targetPositions,blockCopies,
coordinatesCrystal,coordinatesCartesian,coordinatesCrystalEmbedded,newCoordinates,
hostM,hostMinverse,targetPositionsCartesian,M,T,
anchorShift=OptionValue["RotationAnchorShift"],
anchorReference=OptionValue["RotationAnchorReference"],
rotationAxes=OptionValue["RotationAxes"],
rotationMap=OptionValue["RotationMap"],
R,twist
},

(* Checking input *)
blocks=domains/.labelMap;
Scan[InputCheck["CrystalQ",#]&,DeleteDuplicates@blocks];

(* Checking if all blocks have same size *)
blockSizes=($CrystalData[#,"Notes","StructureSize"]/._Missing->{1,1,1})&/@blocks;
If[!SameQ@@blockSizes,Message[SynthesiseStructure::DifferentBlockSizes];Abort[]];
blockSizes=blockSizes[[1]];

(* Checking if output size is valid *)
outputSize={A,B,C};
If[A*B*C!=Length@domains,
Message[SynthesiseStructure::SizeError];Abort[]];
If[(Length@outputSize!=3)||(AnyTrue[outputSize,Positive[#]\[Nand]IntegerQ[#]&]),
Message[SynthesiseStructure::InvalidOutputSize];Abort[]];
If[Total@MapThread[Mod,{outputSize,blockSizes}]=!=0,
Message[SynthesiseStructure::IncompatibleOutputSize];Abort[]];

b=If[TrueQ@OptionValue["UsePlacementBuffer"],2,1];
targetPositions=Flatten[Table[{i,j,k},
{i,0.,b*outputSize[[1]]-1,b*blockSizes[[1]]},
{j,0.,b*outputSize[[2]]-1,b*blockSizes[[2]]},
{k,0.,b*outputSize[[3]]-1,b*blockSizes[[3]]}
],2];

AppendTo[$CrystalData,outputName->$CrystalData@First@blocks];
hostM=GetCrystalMetric[outputName,"ToCartesian"->True];
hostMinverse=Inverse@hostM;
targetPositionsCartesian=hostM.#&/@targetPositions;

If[rotationMap=!=<||>,
R=InputCheck["RotationTransformation",{outputSize,domains},
{anchorShift,anchorReference,rotationMap,rotationAxes},True]
];

blockCopies=$CrystalData/@blocks;
Do[
M=GetCrystalMetric[blocks[[i]],"ToCartesian"->True];
T=TranslationTransform[targetPositions[[i]]];

coordinatesCrystal=blockCopies[[i,"AtomData",All,"FractionalCoordinates"]];
coordinatesCartesian=M.#&/@coordinatesCrystal;
coordinatesCrystalEmbedded=hostMinverse.#&/@coordinatesCartesian;
newCoordinates=T/@coordinatesCrystalEmbedded;

If[rotationMap=!=<||>,
twist=R[domains[[i]],targetPositionsCartesian[[i]]];
coordinatesCartesian=hostM.#&/@newCoordinates;
coordinatesCartesian=twist@coordinatesCartesian;
newCoordinates=hostMinverse.#&/@coordinatesCartesian
];

blockCopies[[i,"AtomData",All,"FractionalCoordinates"]]=newCoordinates,
{i,2,Length@blockCopies}
];

$CrystalData[outputName,"AtomData"]=Flatten@blockCopies[[All,"AtomData"]];
$CrystalData[outputName,"Notes","StructureSize"]=outputSize;

outputName
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
Options@SystematicAbsentQ={
(* Options from 'StructureFactor' *)
"Threshold"->0
};

SyntaxInformation@SystematicAbsentQ={
"ArgumentsPattern"->{_,_,OptionsPattern[]}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
SystematicAbsentQ[group_String,hkl_List,OptionsPattern[]]:=
Block[{sgHM,centring,symop,r,t,\[Delta],crystalQ,\[Lambda],check},

(*---* Checking and processing input *---*)
	If[KeyExistsQ[$CrystalData,group],
	sgHM=$CrystalData[group,"SpaceGroup"],
	sgHM=group];

	sgHM=InputCheck["InterpretSpaceGroup",sgHM];
	InputCheck[hkl,"Integer"];
	
	(* Loading lattice centring vectors and symmetry operations *)
	sgHM=GetSymmetryData[sgHM,"HermannMauguinFull"];
	centring=InputCheck["GetCentringVectors",StringTake[sgHM,1]];
	symop=GetSymmetryOperations[sgHM];

(*---* Procedure for checking systematic absence *---*)
	(* Useful variables *)
	\[Delta]=OptionValue["Threshold"];
	crystalQ=KeyExistsQ[$CrystalData,group];
		(* Ony check structure factor if \[Delta] is numeric *)
		If[NumericQ@\[Delta]&&\[Delta]>0,
		\[Lambda]=$CrystalData[group,"Wavelength"];
		If[MissingQ@\[Lambda],\[Lambda]=1.0],

		crystalQ=False];

check[input_]:=(
(* General absence (centring) *)
	If[!AllTrue[Dot[input,#]&/@centring,IntegerQ],Return@True];

(* Special absence (translation) *)
	r=Equal[input,#]&/@Transpose[Transpose@symop[[All,1]].input];
	t=IntegerQ[input.#]&/@symop[[All,2]];
	If[
	MemberQ[Transpose[{r,t}],{True,False}],
	Return@True];

(* Structure factor below threshold *)
	If[crystalQ,If[Abs[First@
StructureFactor[group,input,\[Lambda],
"Units"->False,
"IgnoreSystematicAbsence"->True]]<\[Delta],
	Return@True]];

(* Not systematically absent reflection *)
	False
	);

(*---* Checking one or more reflections *---*)
If[MatrixQ@hkl,
(* Multiple reflections *)
Reap[Do[Sow[check@hkl[[i]]],
{i,Length@hkl}]][[2,1]],

(* Single reflection input *)
check[hkl]
]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
ToStandardSetting::ext="Invalid extension \[LeftGuillemet]`1`\[RightGuillemet] for this space group.";

SyntaxInformation@ToStandardSetting={
"ArgumentsPattern"->{_,_.}
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
ToStandardSetting[group_String,hkl_List]:=
Block[{operators,equivalents,nonnegative,list},
(* Input check *)
	InputCheck[hkl,"1hkl","Integer"];
	InputCheck["PointSpaceGroupQ",group];

(* Listing all symmetry-equivalents *)
	operators=GetSymmetryOperations@GetLaueClass@group;

	If[Depth@operators==5,
	operators=operators[[All,1]]];

	equivalents=Transpose[#].hkl&/@operators;

(* Selection *)
	nonnegative=Select[equivalents,AllTrue[#,NonNegative]&];

	If[nonnegative!={},
	list=nonnegative,
	list=equivalents];

Last@SortBy[list,{#[[3]]&,#[[2]]&,#[[1]]&}]
]


(* ::Input::Initialization:: *)
ToStandardSetting[input_String,extension_:1]:=
Block[{
sg,fullHM,SG,
mainTargetQ,uniqueInSgQ,mainSetting,
output,temp},

(* Prepartaion *)
output=input;

	(* Check if valid input symbol *)
	InputCheck["PointSpaceGroupQ",input];

	sg=$GroupSymbolRedirect[input];
	fullHM=sg["Name","HermannMauguinFull"];

	temp=Position[$SpaceGroups,fullHM];
	SG=$SpaceGroups[temp[[1,1,1]]];

(* Is target space group main setting? *)
mainTargetQ=Length@First@temp<=3;

(* Is target symbol unique among all space groups? *)
temp=Position[$SpaceGroups,sg["Name","Symbol"]];
If[!Length@DeleteDuplicates@temp[[All,1,1]]===1,
Return@fullHM];

(* Is target symbol unique in its space group? *)
uniqueInSgQ=Length@temp===1;	

	Which[
	(* With cell origin tag? *)
	extension=!=1,
	output=sg["Name","Symbol"]<>":"<>ToString@extension;
	If[!KeyExistsQ[$GroupSymbolRedirect,output],
	Message[ToStandardSetting::ext,extension];Abort[]],

	(* Is target "best candidate" for non-main symbol? *)
	!uniqueInSgQ&&!mainTargetQ,
	mainSetting=First@SG["Setting"];
	If[First@sg["Setting"]===mainSetting,
	output=sg["Name","Symbol"],
	output=fullHM],

	(* Uniquely formatted symbol OR "main symbol"? *)
	uniqueInSgQ||mainTargetQ,
	output=$GroupSymbolRedirect[output]["Name","Symbol"],

	(* If ambiguous, use full Hermman--Mauguin *)
	True,
	output=fullHM
	];

Return@output
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
UnitCellTransformation::command="Transformation command \[LeftGuillemet]`1`\[RightGuillemet] not recognised.";
UnitCellTransformation::invalid="The setting \[LeftGuillemet]`1`\[RightGuillemet] is invalid.";
UnitCellTransformation::invalidSG="The setting \[LeftGuillemet]`1`\[RightGuillemet] is invalid for this particular space group.";
UnitCellTransformation::na="Input not applicable for the `1` system.";
UnitCellTransformation::diffset="Setting mismatch between source, `1`, and target, `2`.";
UnitCellTransformation::sg="Could not interpret \[LeftGuillemet]`1`\[RightGuillemet] as a space group symbol.";
UnitCellTransformation::differentsg="The source space group, \[LeftGuillemet]`1`\[RightGuillemet] (no.`2`), and the target space group, \[LeftGuillemet]`3`\[RightGuillemet] (no. `4`), are not the same.";
UnitCellTransformation::support="Transformation of the input is not supported.";
UnitCellTransformation::target="Could not determine target space group.";
UnitCellTransformation::one="This space group, \[LeftGuillemet]`1`\[RightGuillemet] (no. `2`), has no alternative representations.";

Options@UnitCellTransformation={
"ReturnP"->False,
"CustomP"->False
};


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
UnitCellTransformation[crystal_String,userinput___]:=
Block[{
(*---* 1. Input check and preparation *---*)
(* 1.A. Load crystal metric, space group and crystal system *)
G,G0,sg,sg0,sourceSG,fullHM,posSG,SG,fullHMs,sourceSGnumber,system,mainSourceQ,
P,P1,P2,p,
centList,axpList,tetList,hex3List,hexList,monoList,
sourceSetting,targetSetting,
(* 1.B. Process input syntax and options *)
inputRules,inputString,returnP,customP,
(* 1.C. Interpret space group from input *)
targetSG,needtargetSG,
(* 1.D Process source setting *)
sourceCentring,sourceCell,notes,relevantNotes,sourceO,
(* 1.E Process target setting *)
allowed,cmds,na,
(* 1.F. Validating input values *)
targetCentring,targetAxis,targetCC,targetAP,targetCell,targetRS,
targetO,
(* 1.G. Check setting constraints on certain space groups *)
(* 1.H. Determining new space group symbol *)
targetSGnumber,
(* 1.I. Common transformation procedures *)
procedureCellCentring,procedureCellOrigin,shift,
(*---* 2. Determining correct transformation matrices *---*)
(* 2.A. Triclinic *)
(* 2.B. Monoclinic *)
cmd,Q,
(* 2.C. Orthorhombic *)
(* 2.D. Tetragonal *)
(* 2.E. Hexagonal crystal family *)
sourceRS,M,
(* 2.F. Cubic *)
(*---* 3. Carrying out transformations *---*)
(* 3.A. Preparations *)
q,newlattice,
(* 3.B. Transforming coordinates and ADPs *)
xyz,newxyz,adps,U,n0,n,newU,u,
(*---* 4. Overwriting entry in $CrystalData *---*)
targetFullHM,
(*---* 5. Display *---*)

(* Dummy variables *)
temp,x,y,i},

(*---* 1. Input check and preparation *---*)
(* 1.A. Load crystal metric, space group and crystal system *)
	InputCheck["GetCrystalSpaceGroup",crystal];
	G=G0=GetCrystalMetric@crystal;
	sg=sg0=$CrystalData[crystal,"SpaceGroup"];
	sg=$GroupSymbolRedirect[sg];
	sourceSG=fullHM=sg["Name","HermannMauguinFull"];
	posSG=Position[$SpaceGroups,fullHM];
	SG=$SpaceGroups[posSG[[1,1,1]]];
	fullHMs=SG["Name","HermannMauguinFull"];
	sourceSGnumber=SG["SpaceGroupNumber"];
	system=SG["CrystalSystem"];
	mainSourceQ=Length@First@posSG<=3;

	(* Default transformation matrices and origin shift *)
	P=P1=P2=IdentityMatrix[3];
	p={0,0,0};

	(* Frequently used lists *)
	centList={
	"P","A","B","C","F","I"};
	axpList={
	"abc","ba-c","cab","-cba","bca","a-cb"};
	tetList={
	"P","I","C1","F1","C2","F2"};
	hex3List={
	"R1","R2","R3"};
	hexList={
	"C1","C2","C3","H1","H2","H3","D1","D2"};
	monoList={
	"Cb1","Cb2","Cb3","Ac1","Ac2","Ac3"};

	(* Miscellaneous *)
	sourceSetting=targetSetting=<||>;

(* 1.B. Process input syntax and options *)
	inputRules=Association@Cases[{userinput},_Rule];
	inputString=Cases[{userinput},_String];
		(* Check *)
		Which[
		#===0,Null,
		#===1,inputString=First@inputString,
		True,
			Message[UnitCellTransformation::invalid,inputString];
			Abort[]
		]&@Length@inputString;

	(* Save option settings *)
	returnP=inputRules["ReturnP"];
	customP=Lookup[inputRules,"CustomP",{}];

	(* Remove options from other settings *)
	KeyDropFrom[inputRules,
	{"ReturnP","CustomP","CustomSymbol"}];

(* 1.C. Interpret space group from input *)
	(* i. No input commands -- prompt dialogue/UI (TODO) *)
	
	(* ii. Custom transformation matrix *)
	If[customP=!={},P1=customP;
	If[!MatrixQ@P1,(* Check *)
	Message[UnitCellTransformation::invalid,P1];Abort[]];
	targetSG=sourceSG(* Use same space group *);
	Goto["MetricTransformation"]];

	(* iii. Interpret string as target space group *)
	If[inputString=!={},

		(* Check if valid space group symbol *)
		If[MissingQ@$GroupSymbolRedirect[inputString],
		Message[UnitCellTransformation::sg,inputString];
		Abort[]];

		targetSG=inputString];

	(* iv. Parse input as setting commands *)
	If[inputRules=!=<||>,
	needtargetSG=True,
		(* Check whether Hall symbol has been used *)
		temp=Position[$SpaceGroups,targetSG];
		If[temp!={},
		If[temp[[1,-1,1]]==="HallString",
		targetSG=First[$SpaceGroups@@@
		Most@First@temp]["Name","Symbol"]]]
	];

(* 1.D Process source setting *)
	(* i. Load source setting from source space group *)
	sourceSetting=$GroupSymbolRedirect[sourceSG]["Setting"];

	(* ii. Check if space group has alternative settings *)
	If[sourceSetting===<||>&&
	(* Exception: Special multiple cells *)
	!MemberQ[{"Tetragonal","Hexagonal"},system],
		Message[UnitCellTransformation::one,
		GetSymmetryData[sourceSG,"Symbol"],
		sourceSGnumber];
		Abort[]];

	(* iii. Check cell origin from symbol *)
	temp=StringTake[sg0,-2];
	Which[
	temp===":2",AppendTo[sourceSetting,"CellOrigin"->2]
	];

	(* iv. Checking 'Notes' for info on input setting *)
	notes=Lookup[$CrystalData@crystal,"Notes",<||>];
	relevantNotes=notes[[{
	"RhombohedralSetting","MultipleCell","CellCentring","CellOrigin"}]];
	{sourceRS,sourceCell,sourceCentring,sourceO}=Values@relevantNotes;
	AppendTo[sourceSetting,DeleteMissing@relevantNotes];

(* 1.E Process target setting *)
	If[inputRules=!=<||>,
	(* a. Setting commands given in association *)
		targetSetting=inputRules;

			(* Checking commands *)
			allowed={
			"UniqueAxis","CellChoice","AxisPermutation",
			"CellCentring","MultipleCell",
			"RhombohedralSetting","CellOrigin"};
			cmds=Keys@targetSetting;
			na=Complement[cmds,allowed];
				If[na!={},
				Message[UnitCellTransformation::command,First@na];
				Abort[]];

			(* Checking usefullness of commands *)
			temp=Which[
			system==="Triclinic",
				{"CellCentring"},
			system==="Monoclinic",
				{"UniqueAxis","CellChoice"},
			system==="Orthorhombic",
				{"AxisPermutation","CellCentring","CellOrigin"},
			system==="Tetragonal",
				{"CellCentring","MultipleCell","CellOrigin"},
			system==="Trigonal",
				{"MultipleCell","RhombohedralSetting"},
			system==="Hexagonal",
				{"MultipleCell"},
			system==="Cubic",
				{"CellCentring","CellOrigin"}
			];

			targetSetting=DeleteMissing@targetSetting[[temp]];
			If[targetSetting===<||>,
			Message[UnitCellTransformation::na,
			ToLowerCase@system];
			Abort[]],

		(* b. Setting extracted from target space group *)
		If[!ValueQ@targetSG,targetSG=sourceSG];
		targetSetting=$GroupSymbolRedirect[targetSG]["Setting"]
	];

	(* Supply with current settings if unspecified in input *)
	targetSetting=Append[sourceSetting,targetSetting];
		(* Exception: No 'RhombohedralSetting' for 'R' cell *)
		If[targetSetting["MultipleCell"]==="R",
		KeyDropFrom[targetSetting,"RhombohedralSetting"]]; 

(* 1.F. Validating input values *)
	(* i. 'CellCentring' *)
	targetCentring=targetSetting["CellCentring"];
	If[!MissingQ@targetCentring,
	If[!MemberQ[centList,
	targetCentring],
	Message[UnitCellTransformation::invalid,
	targetCentring];Abort[]
	]];

	(* ii. 'UniqueAxis' *)
	targetAxis=targetSetting["UniqueAxis"];
	If[!MissingQ@targetAxis,
	If[!MemberQ[{"a","b","c"},targetAxis],
	Message[UnitCellTransformation::invalid,
	targetAxis];Abort[]];
	targetAxis=ToLowerCase@targetAxis
	];

	(* iii. 'CellChoice' *)
	targetCC=targetSetting["CellChoice"];
	If[!MissingQ@targetCC,
	If[!MemberQ[{1,2,3},targetCC],
	Message[UnitCellTransformation::invalid,
	targetCC];Abort[]]
	];

	(* iv. 'AxisPermutation' *)
	targetAP=targetSetting["AxisPermutation"];
	If[!MissingQ@targetAP,
	(* Check if string *)
	If[!StringQ@targetAP,
	Message[UnitCellTransformation::invalid,
	targetAP];Abort[]];
	(* Support for various input forms *)
	targetAP=ToLowerCase@StringReplace[
	targetAP,{"\!\(\*OverscriptBox[\(c\), \(_\)]\)"->"-c","\!\(\*OverscriptBox[\(C\), \(_\)]\)"->"-C"}];
	(* Check setting value *)
	If[!MemberQ[axpList,targetAP],
	Message[UnitCellTransformation::invalid,
	targetSetting["AxisPermutation"]];Abort[]];
	(* Update 'targetSetting' *)
	targetSetting["AxisPermutation"]=targetAP];

	(* v. 'MultipleCell' *)
	targetCell=targetSetting["MultipleCell"];
	If[!MissingQ@targetCell,
	(* Depending on crystal system... *)
	temp=Which[
	system==="Tetragonal",tetList,
	system==="Trigonal"||system==="Hexagonal",
		Flatten@Join[
		{"P","R"},hex3List,hexList,monoList],
	True,
		Message[UnitCellTransformation::na,
		ToLowerCase@system];Abort[]];
	(* Check setting value *)
	If[!MemberQ[temp,targetCell],
	Message[UnitCellTransformation::invalidSG,
	targetCell];Abort[]]
	];

	(* vi. 'RhombohedralSetting' *)
	targetRS=inputRules["RhombohedralSetting"];
	If[!MissingQ@targetRS,
	targetRS=ToLowerCase@targetRS;
	(* Check setting value *)
	If[!MemberQ[{"obverse","reverse"},targetRS],
	Message[UnitCellTransformation::invalid,
	targetSetting["RhombohedralSetting"]];Abort[]];
	(* Update 'targetSetting' *)
	targetSetting["RhombohedralSetting"]=targetRS
	];

	(* vii. 'CellOrigin' *)
	targetO=targetSetting["CellOrigin"];
	If[!MissingQ@targetO,
	If[!MemberQ[{1,2},targetO],
	Message[UnitCellTransformation::invalid,
	targetO];Abort[]]
	];

(* 1.G. Check setting constraints on certain space groups *)
	(* i. Target setting must be a subset of space group settings *)
	If[!(SubsetQ@@Keys/@{SG["Setting"],
Which[(* Exceptions: Cell centring and special multiple cells *)
MemberQ[{"Tetragonal","Hexagonal"},system],
	KeyDrop[targetSetting,"MultipleCell"],
MemberQ[{"Monoclinic","Orthorhombic","Tetragonal","Cubic"},system],
	KeyDrop[targetSetting,"CellCentring"],
True,
	targetSetting]}
	),Message[UnitCellTransformation::diffset,
		Keys@sourceSetting,Keys@targetSetting];Abort[]
	];

(* 1.H. Determining new space group symbol *)
	(* i. Determine the target space group from setting if needed *)
		If[needtargetSG,
		(* Main entry? *)
		If[Sort@SG["Setting"]===Sort@targetSetting,
		targetSG=SG["Name","HermannMauguinFull"],
			(* Special multiple cell? *)
			If[MemberQ[Flatten@Join[
			{tetList,{"R"},hex3List,hexList,monoList}],
			targetCell],
			targetSG=sourceSG;
			targetSGnumber=sourceSGnumber;
			Goto["CheckSGformat"]
			];

		(* Check alternatives *)
		temp=SG[["AlternativeSettings"]];
		i=1;
		While[True,
		If[i>Length@temp,
		Message[UnitCellTransformation::target];Abort[]];
		If[Sort@temp[[i,"Setting"]]===Sort@targetSetting,
		targetSG=temp[[i,"Name","HermannMauguinFull"]];
		Break[]];
		i++]
		]];

	(* ii. Check if source and target are same space group *)
	targetSGnumber=GetSymmetryData[targetSG,"SpaceGroupNumber"];

	If[sourceSGnumber!=targetSGnumber,
	Message[UnitCellTransformation::differentsg,
	GetSymmetryData[sourceSG,"Symbol"],sourceSGnumber,
	GetSymmetryData[targetSG,"Symbol"],targetSGnumber];
	Abort[]];

	(* iii. Check whether to use formatted symbol *)
	Label["CheckSGformat"];
	targetSG=ToStandardSetting@targetSG;

(* 1.I. Common transformation procedures *)
	(* 'CellCentring' *)
	procedureCellCentring:=(
	If[!MissingQ@targetCentring,
	If[MissingQ@sourceCentring||!ValueQ@sourceCentring,
	sourceCentring=StringTake[fullHM,1]];

	Which[
	(* No change *)
	sourceCentring===targetCentring,
	Null,

	(* Transformation from 'P' *)
	sourceCentring==="P",
	P1=Inverse@$TransformationMatrices[
	targetCentring<>"_to_P"],

	(* Transformation to 'P' *)
	targetCentring==="P",
	P1=$TransformationMatrices[sourceCentring<>"_to_P"],

	(* Transformation via 'P' *)
	True,
	P1=$TransformationMatrices[sourceCentring<>"_to_P"];
	P2=Inverse@$TransformationMatrices[targetCentring<>"_to_P"]
	]
	]);

	(* 'CellOrigin' *)
	procedureCellOrigin:=(
	If[KeyExistsQ[targetSetting,"CellOrigin"],
	temp={sourceSetting["CellOrigin"],targetO};
	shift=SG[["AlternativeSettings",1,"OriginShift"]];
	If[!MissingQ@targetO,
	(* Get space group origin shift (shift vector) *)
	Which[
	temp==={1,1},Null,
	temp==={2,2},Null,
	temp==={1,2},p=shift,
	temp==={2,1},p=-shift]
	]
	]);


(*---* 2. Determining correct transformation matrices *---*)
Goto[system];

(*-- 2.A. Triclinic --*)
Label["Triclinic"];

	(* i. 'CellCentring' *)
	procedureCellCentring;

	(* Preparations done *)
	Goto["MetricTransformation"];


(*-- 2.B. Monoclinic --*)
Label["Monoclinic"];

	(* Prepartations *)
		(* If target axis not given, use same axis as source *)
		If[MissingQ@targetAxis||!ValueQ@targetAxis,
		targetAxis=sourceSetting["UniqueAxis"]];

	(* i. 'UniqueAxis' *)
		(* Target unique axis transformation *)
		allowed={
		"UniqueAxisB_to_C",
		"UniqueAxisB_to_A",
		"UniqueAxisC_to_A"};

		temp=ToUpperCase/@{sourceSetting["UniqueAxis"],targetAxis};
	
		(* Check whether P or Q (inverse) is needed *)
		cmd="UniqueAxis"<>temp[[1]]<>"_to_"<>temp[[2]];
		If[MemberQ[allowed,cmd],
		Q=False,
		Q=True;
			cmd="UniqueAxis"<>temp[[2]]<>"_to_"<>temp[[1]]
		];

		(* Target unique axis transformation *)
		If[!SameQ@@temp,
		P1=$TransformationMatrices[cmd];
		If[Q,P1=Inverse@P1]
		];

	(* ii. 'CellChoice' *)
		(* Check if cell choice is an available setting *)		
		If[!KeyExistsQ[SG["Setting"],"CellChoice"]
		&&KeyExistsQ[targetSetting,"CellChoice"],
		Message[UnitCellTransformation::command,"CellChoice"];
		Abort[],

		If[!MissingQ@targetCC,
		(* Matrix for checking transformation signature *)
		temp={{0,1,-1},{-1,0,1},{1,-1,0}};
		temp=temp[[sg["Setting","CellChoice"],targetCC]];
		(* Use regular, none, or inverse transformation *)
		Which[
		temp===-1,Q=True,
		temp===0,Goto["MetricTransformation"],
		temp===1,Q=False
		];
		cmd="UniqueAxis"<>ToUpperCase@targetAxis
		<>"_CellChoice+1";
		P2=$TransformationMatrices[cmd];
		If[Q,P2=Inverse@P2];
		]];

	(* Preparations done *)
	Goto["MetricTransformation"];


(*-- 2.C. Orthorhombic --*)
Label["Orthorhombic"];

	(* i. 'AxisPermutation' *)
		If[!MissingQ@targetAP,
		(* Check which matrices are needed *)
		temp=ToUpperCase/@
		{sourceSetting["AxisPermutation"],targetAP};
		If[SameQ@@temp,
		(* Same axis permutation -- no transform required *)
		Null,
		If[!MemberQ[temp,"DEF"],
		(* Chain through 'abc' *)
		x=temp[[1]]<>"_to_DEF";
		y=temp[[2]]<>"_to_DEF";
		P1=$TransformationMatrices[x];
		P2=Inverse@$TransformationMatrices[y],
		(* OR use one of the six transformations *)
		cmd=temp[[1]]<>"_to_"<>temp[[2]];
		If[KeyExistsQ[$TransformationMatrices,cmd],
		P1=$TransformationMatrices[cmd],
		P1=Inverse@$TransformationMatrices[
	temp[[2]]<>"_to_"<>temp[[1]]]
		]
		]
		]
		];

	(* ii. 'CellOrigin' *)
	procedureCellOrigin;

	(* Preparations done *)
	Goto["MetricTransformation"];


(*-- 2.D. Tetragonal --*)
Label["Tetragonal"];

	(* Checks and updates *)
	If[MissingQ@sourceCentring||!ValueQ@sourceCentring,
	sourceCentring=StringTake[
	sg["Name","HermannMauguinShort"],1]];

	If[MissingQ@targetCentring||!ValueQ@targetCentring,
	targetCentring=StringTake[
	SG["Name","HermannMauguinShort"],1]];

	(* i. 'CellCentring' *)
	procedureCellCentring;

	(* ii. 'MultipleCell' *)
	If[KeyExistsQ[targetSetting,"MultipleCell"],
	(* If already transformed, transform to 'P' or 'I' *)
	If[MemberQ[tetList,sourceCell],
	P1=Inverse@$TransformationMatrices[
	"TetragonalProjection"<>StringTake[sourceCell,-1]]
	];

	(* Transforming to target cell *)
	If[!MemberQ[{"P","I"},targetCell],
	P2=$TransformationMatrices[
	"TetragonalProjection"<>StringTake[targetCell,-1]]
	]
	];
	
	(* iii. 'CellOrigin' *)
	procedureCellOrigin;

	(* Preparations done *)
	Goto["MetricTransformation"];


(*-- 2.E. Hexagonal crystal family --*)
Label["Trigonal"];
Label["Hexagonal"];

	(* Checks and updates *)
		(* Reinstate 'MultipleCell' label *)
		If[MissingQ@sourceCell||!ValueQ@sourceCell,
		sourceCell=sourceSetting["MultipleCell"]];

		(* Check rhombohedral source setting *)
		If[MissingQ@sourceRS||!ValueQ@sourceRS,
		(* Assume main space group setting *)
		sourceRS="obverse"];

		(* Check rhombohedral target setting *)
		If[MissingQ@targetRS&&MemberQ[hex3List,targetCell],targetRS="obverse"];
		If[!MissingQ@targetRS&&targetCell==="R",targetCell="R1"];

	Which[
	(* A. Rhombohedral space group *)
	SubsetQ[{146,148,155,160,161,166,167},
	{sourceSGnumber,targetSGnumber}],

		(* Check target command *)
		If[!MemberQ[Flatten@Join[
	{{"R"},hex3List,monoList}],targetCell],
		Message[UnitCellTransformation::invalidSG,
			targetCell];Abort[]];

		Which[
		(* a. Rhombohedral source *)
		sourceCell==="R",
		Which[
		(* a.1. Rhombohedral target *)
		targetCell==="R"&&MissingQ@targetRS,
			Goto["MatricesDone"],

		(* a.2. Triple hexagonal target *)
		MemberQ[hex3List,targetCell],
		P1=$TransformationMatrices[
		"Rhombohedral_to_"<>targetCell<>"_"<>targetRS],

		(* a.3. Monoclinic target *)
		True,
		P1=$TransformationMatrices[
		"Rhombohedral_to_Monoclinic_"<>targetCell]
		],

		(* b. Regular hexagonal (R1, R2, R3) source *)
		MemberQ[hex3List,sourceCell],
		Which[
		(* b.1. Rhombohedral target *)
		targetCell==="R",
		P1=Inverse@$TransformationMatrices[
		"Rhombohedral_to_"<>sourceCell<>"_"<>sourceRS],

		(* b.2. Triple hexagonal target *)
		MemberQ[hex3List,targetCell],
		If[sourceCell===targetCell&&sourceRS===targetRS,
		Goto["MatricesDone"]];
		P1=Inverse@$TransformationMatrices[
		"Rhombohedral_to_"<>sourceCell<>"_"<>sourceRS];
		P2=$TransformationMatrices[
		"Rhombohedral_to_"<>targetCell<>"_"<>targetRS],

		(* b.3. Monoclinic target *)
		True,
		If[sourceCell==="R1"&&sourceRS==="obverse",
		P1=$TransformationMatrices[
		"TripleHexagonal_to_Monoclinic_"<>targetCell],
		P1=Inverse@$TransformationMatrices[
		"Rhombohedral_to_"sourceCell<>"_"<>sourceRS];
		P2=$TransformationMatrices[
		"Rhombohedral_to_Monoclinic_"<>targetCell]]
		],


		(* c. Monoclinic source *)
		MemberQ[monoList,sourceCell],
		Which[
		(* c.1. Rhombohedral target *)
		targetCell==="R",
		P1=Inverse@$TransformationMatrices[
		"Rhombohedral_to_Monoclinic_"<>sourceCell],

		(* c.2. Regular hexagonal (R1, R2, R3)l target *)
		MemberQ[hex3List,targetCell],
		P1=Inverse@$TransformationMatrices[
		"Rhombohedral_to_Monoclinic_"<>sourceCell];
		P2=$TransformationMatrices[
		"Rhombohedral_to_"<>targetCell<>"_"<>targetRS],

		(* c.3. Monoclinic target *)
		True,
		If[sourceCell===targetCell,Goto["MatricesDone"]];
		P1=Inverse@$TransformationMatrices[
		"Rhombohedral_to_Monoclinic_"<>sourceCell];
		P2=$TransformationMatrices[
		"Rhombohedral_to_Monoclinic_"<>targetCell]
		]
		],
	

	(* B. Transformation of the hexagonal lattice *)
	True,
	
	(* Check if target is valid *)
	If[!MemberQ[Flatten@Join[
	{{"P"},hexList,monoList}],targetCell],
		Message[UnitCellTransformation::invalidSG,
			targetCell];Abort[]];

	M=Flatten@DeleteCases[
	StringCases[Keys@$TransformationMatrices,
	StartOfString~~"Hexagonal"~~__],{}];

	If[MissingQ@sourceCell,sourceCell="P"];

		Which[
		(* a. Primitive hexagonal source *)
		sourceCell==="P",
		Which[
		(* a.1. No change *)
		targetCell==="P",Goto["MatricesDone"],
		(* a.2. Special hexagonal target *)
		MemberQ[hexList,targetCell],
		temp=Select[M,StringContainsQ[#,targetCell]&];
		P1=$TransformationMatrices@@temp
		],

		(* b. Special hexagonal source *)
		MemberQ[hexList,sourceCell],
		Which[
		(* b.1. To primitive hexagonal cell *)
		targetCell==="P",
		temp=Select[M,StringContainsQ[#,targetCell]&];
		P1=Inverse@$TransformationMatrices@@temp,

		(* b.2. To special hexagonal target *)
		MemberQ[hexList,targetCell],
		If[sourceCell===targetCell,Goto["MatricesDone"]];
		temp=Select[M,StringContainsQ[#,sourceCell]&];
		P1=Inverse[$TransformationMatrices@@temp];
		temp=Select[M,StringContainsQ[#,targetCell]&];
		P2=$TransformationMatrices@@temp
		]
		]
	];
	Label["MatricesDone"];

	(* Preparations done *)
	Goto["MetricTransformation"];


(*-- 2.F. Cubic --*)
Label["Cubic"];

	(* i. 'CellCentring' *) 
	procedureCellCentring;

	(* ii. 'CellOrigin' *)
	procedureCellOrigin;

	(* Preparations done *)
	Goto["MetricTransformation"];


(*---* 3. Carrying out transformations *---*)
(* 3.A. Preparations *)
Label["MetricTransformation"];
	P=P1.P2;
		If[returnP===All,Return@{P1,P2}];
		If[returnP,Return@P];
	Q=Inverse@P;
	q=-Q.p;
	G=Transpose[P].G.P;

	newlattice=Association@Thread[
{"a","b","c","\[Alpha]","\[Beta]","\[Gamma]"}
->GetLatticeParameters[G,"Units"->True]];

(* 3.B. Transforming coordinates and ADPs *)
	(* Fractional coordinates *)
	xyz=$CrystalData[[crystal,
	"AtomData",All,"FractionalCoordinates"]];
	newxyz=Chop[FractionalPart[Dot[Q,#]+q]&/@xyz];

	(* Atomic displacement parameters *)
	adps=Lookup[$CrystalData[crystal,"AtomData"],
	"DisplacementParameters",0.];
	U={{#1,#4,#5},{#4,#2,#6},{#5,#6,#3}}&@@#&/@adps;

	(* Preparing diagonal 'N' matrices *)
	(* References:
	https://doi.org/10.1107/S0108767311018216
	https://doi.org/10.1107/S0021889802008580 *)
	n0=DiagonalMatrix@Sqrt@Diagonal@Inverse@G0;
	n=Inverse@DiagonalMatrix@Sqrt@Diagonal@Inverse@G;

	(* Transforming ADPs *)
	newU={};
	Do[
	u=U[[i]];
	If[MatrixQ[u],
	temp=n0.u.Transpose[n0];
	temp=Q.temp.Transpose[Q];
	temp=Chop[n.temp.Transpose[n]];
	temp=Part[temp,#/.List->Sequence]&/@
	{{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}},
	temp=u];
	AppendTo[newU,temp],
	{i,Length@adps}
	];


(*---* 4. Overwriting entry in $CrystalData *---*)
	(* Determine new space group symbol automatically *)
		(* Remove old cell tag *)
		If[StringTake[targetSG,{-2}]===":",
		targetSG=StringTake[targetSG,{1,-3}]];

		(* Update space group symbol *)
		If[(targetO===2||targetSetting["CellOrigin"]===2)&&
		!StringContainsQ[targetSG,":2"],
		targetSG=targetSG<>":2"];

		(* Get full Hermann\[Dash]Mauguin symbol *)
		targetFullHM=GetSymmetryData[targetSG,"HermannMauguinFull"];

		(* Set new symbol *)
		$CrystalData[crystal,"SpaceGroup"]=targetSG;

	(* New lattice parameters *)
	AppendTo[$CrystalData[crystal,
	"LatticeParameters"],newlattice];

	(* New fractional coordinates *)
	$CrystalData[[crystal,"AtomData",All,
	"FractionalCoordinates"]]=newxyz;

	(* New ADPs *)
	If[!AllTrue[N@newU,#==0.&],
	$CrystalData[[crystal,"AtomData",All,
	"DisplacementParameters"]]=newU];

	(* Updating 'Notes' and adding space group notes if needed *)
	If[
	(* a. Exception: default space group setting *)
	Sort@Values@SG["Setting"]===
	Sort@DeleteMissing@{targetCell,targetRS,targetCentring},
	KeyDropFrom[notes,{"MultipleCell","RhombohedralSetting","CellCentring"}],


	(* b. Write both multiple cell and rhombohedral setting *)
	If[KeyExistsQ[targetSetting,"MultipleCell"]&&
	!MemberQ[{"P","I"},targetCell],
	AppendTo[notes,"MultipleCell"->targetCell]];
	
	If[MemberQ[hex3List,targetCell],
	AppendTo[notes,"RhombohedralSetting"->targetRS],
	KeyDropFrom[notes,"RhombohedralSetting"]];

	If[!MissingQ@targetCentring&&ValueQ@targetCentring&&
	(* No need if centring is in space group symbol *)
	StringTake[targetFullHM,1]!=targetCentring,
	AppendTo[notes,"CellCentring"->targetCentring]]
	];

	(* Writing over 'Notes' or removing if empty *)
	If[notes=!=<||>,
	$CrystalData[crystal,"Notes"]=Sort@notes,
	KeyDropFrom[$CrystalData[crystal],"Notes"]
	];


(*---* 5. Display *---*)
Label["End"];

KeyValueMap[
If[#1==="AtomData",
"AtomData"->Shallow[#2,1],
#1->#2]&,
$CrystalData[crystal]]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
$CrystalData::type="The atomic displacement type \[LeftGuillemet]`1`\[RightGuillemet] is not recognised.";
$CrystalData::missing="Could not find \[LeftGuillemet]`1`\[RightGuillemet].";


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
$CrystalData:=$CrystalData=Import[
FileNameJoin[{$MaXrdPath,"UserData","CrystalData.m"}],"Package"];


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
$LaueClasses::notLaue="\[LeftGuillemet]`1`\[RightGuillemet] is not a valid Laue class.";


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
$LaueClasses:=$LaueClasses=$PointGroups[[{"-1","2/m","mmm","4/m","4/mmm","-3","-3m","6/m","6/mmm","m-3","m-3m"}]];


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
$PeriodicTable::data="No data found on \[LeftGuillemet]`1`\[RightGuillemet].";


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
$PeriodicTable:=$PeriodicTable=Import[FileNameJoin[{$MaXrdPath,"Core","Data","PeriodicTable.m"}],"Package"];


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
$PointGroups::symbol="No data found on \[LeftGuillemet]`1`\[RightGuillemet].";


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
$PointGroups:=$PointGroups=Import[FileNameJoin[{$MaXrdPath,"Core","Data","PointGroups.m"}],"Package"];


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
$SpaceGroups::symbol="No data found on \[LeftGuillemet]`1`\[RightGuillemet].";


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
$SpaceGroups:=$SpaceGroups=Import[FileNameJoin[{$MaXrdPath,"Core","Data","SpaceGroups.m"}],"Package"];


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
$TransformationMatrices::incompatible="Incompatible transformation request.";


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
$TransformationMatrices:=$TransformationMatrices=Import[FileNameJoin[{$MaXrdPath,"Core","Data","TransformationMatrices.m"}],"Package"];


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
$MaXrdChangelog::fileMissing="Cannot find the `1` file (`2`).";


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
$MaXrdChangelog:=Block[{
dir=$MaXrdPath,
pacletFile,packletVersion,
packageSymbols,
changelogFile,log,current,content,new,t,title,post,
temp},
(* Load current version -- could be undeployed *)
	pacletFile=FileNames["PacletInfo.m",{#}]
	&/@{dir,ParentDirectory@dir};
	pacletFile=Cases[pacletFile,_String,2,1];

	If[pacletFile==={},
	Message[$MaXrdChangelog::fileMissing,"paclet",pacletFile];
	Abort[],
	pacletFile=pacletFile[[1]]
	];

	packletVersion=(Association@@Import@pacletFile)[Version];

(* Latest entry in the changelog *)
	changelogFile=FileNames["Changelog.md",{dir},2][[1]];
	If[!FileExistsQ@changelogFile,
	Message[$MaXrdChangelog::fileMissing,"changelog",changelogFile];
	Abort[]];
	log=Import[changelogFile,"Text"];
	
	current=Check[StringCases[log,Shortest[
"## Version "~~packletVersion~~"\n"~~news__~~
{"\n\n\n##",EndOfString}]:>news],
	Abort[]];

	(* Content *)
	content=StringSplit[First@current,"\n"..];
	content=StringDrop[#,2]&/@content;
	content=TextCell[Row[{#}],"Item"]&/@content;

	(* Format function names *)
	packageSymbols=First/@Cases[
	Flatten@First@$MaXrdFunctions,_Hyperlink];
	new={};
	Do[
	temp=content[[i,1,1,1]];
	t={StartOfString,"",".",",","!"};
	AppendTo[new,
	Row@StringSplit[temp,{
	Shortest[a:t~~"**"~~x__~~"**"~~b:t]:>
	Style[a<>x<>b,Bold],
	Shortest[a:t~~"_"~~x__~~"_"~~b:t]:>
	Style[a<>x<>b,Italic],
	Shortest[a:t~~"*"~~x__~~"*"~~b:t]:>
	Style[a<>x<>b,Italic],
	Shortest[a:t~~"`"~~x__~~"`"~~b:t]:>
	If[MemberQ[packageSymbols,x],
	Hyperlink[x,"paclet:/MaXrd/Ref/"<>x],
	Style[a<>x<>b,"Program"]]
	}]],
	{i,Length@content}];

	new=(
first=#[[1,1]];
If[StringQ@first&&StringTake[first,2]==="# ",
TextCell[
ReplacePart[#,{1,1}->StringDrop[first,2]],
"Subsubsection"],
TextCell[#,"Item"]
])&/@new;

(* Title *)
	title=TextCell@Row[Style[#,"Text",
FontSize->24,
FontFamily->"Source Sans Pro",
Italic]&/@{
"Changelog for the MaXrd",
" package \[LongDash] version "<>packletVersion}
];

(* Log file link *)
	MaXrd`Private`$changelogFile=changelogFile;
	post=TextCell[
	Row[{
	Mouseover[#1,#2]&@@(Button["Click here",
SystemOpen@MaXrd`Private`$changelogFile,
Appearance->"Frameless",
BaseStyle->{"GenericButton",#}]&/@{
RGBColor[0.269993,0.308507,0.6],
RGBColor[0.823529411764706,0.490196078431373,0.133333333333333]}),
" to see the complete changelog of all versions."
	}],
	"Item"
	];

(* Print all cells *)
CreateDocument[CellGroup@Join[{title},new,{Row[{}]},{post}],CellGrouping->Manual]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
$MaXrdFunctions:=$MaXrdFunctions=
Block[{subContexts,packagefunctions,packagef,hyperlinks},
subContexts=Select[Contexts["MaXrd`*"],!StringContainsQ[#,"Private"]&];
packagefunctions=#->Names[#<>"*"]&/@subContexts;
packagef=Sort@Flatten@packagefunctions[[All,2]];
hyperlinks=Hyperlink[#,"paclet:/MaXrd/Ref/"<>#]&/@packagef;
Multicolumn[hyperlinks,2]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
$MaXrdPath::MissingDefinitionsFile="Unable to locate the package definition file.";


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
$MaXrdPath:=$MaXrdPath=Block[{files,prioritisedFile,definitionFile},
files=FileNames["MaXrd/Core/Definitions.m",$Path];
If[files==={},
Message[$MaXrdPath::MissingDefinitionsFile];Abort[]];
prioritisedFile=Select[files,StringContainsQ[#,FileNameJoin[
{"Mathematica","Applications","MaXrd","Core","Definitions.m"}]]&];
definitionFile=If[prioritisedFile=!={},
prioritisedFile,files][[1]];
DirectoryName[definitionFile,2]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
(* Messages, options, attributes and syntax information *)


(* ::Input::Initialization:: *)
$GroupSymbolRedirect::missing="Could not find \[LeftGuillemet]`1`\[RightGuillemet].";


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
$GroupSymbolRedirect:=$GroupSymbolRedirect=Import[FileNameJoin[{$MaXrdPath,"Core","Data","GroupSymbolRedirect.m"}],"Package"];


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
$MaXrdVersion:=Block[{dir,p},
(* Load current version *)
	dir=$MaXrdPath;
	p=FileNameJoin[{dir,"PacletInfo.m"}];
	First@StringCases[Import[p,"String"],
	Shortest[Whitespace~~"Version -> \""~~v__~~"\""]:>v]
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
EndPackage[];
